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A  Quasi-Equilibrium  Kinetic  Model  of  HF  Air  Breakdown 


I.  INTRODUCTION. 

A.  BACKGROUND. 

This  project  was  initiated  with  the  express  purpose  of  developing  a  model  of  air 
breakdown  in  apertures  subjected  to  high-frequency  (HF)  electromagnetic  pulses. 
High— altitude  electromagnetic  pulse  (HEMP)  threat  levels  and  frequency  content  are 
determined  primarily  by  the  device  rise  time,  and  for  nominal  parameters,  are  in 
the  regime  of  30  kV/m  with  significant  energy  content  up  to  lO’s  of  MHz 
[Ref.  1].  The  requirement  for  the  project  was  inspired  by  the  observation  of 
ubiquitous  shot— to— shot  variations  in  the  response  of  shielded  systems  to  EMP  in 
simulation  experiments,  leading  to  the  hypothesis  that  air  breakdown  was  occurring 
in  the  seams  of  apertures.  Methods  to  solve  the  field  distributions  near  an 
aperture  are  well  known,  at  least  in  the  linear  regime  before  the  air  breakdown. 
The  challenge  is  to  correctly  model  the  highly  nonlinear  breakdown  processes  for 
an  EMP— like  driving  field  pulse.  Because  of  the  nonlinearity  of  the  reactions,  a 
time-domain  approach  is  needed.  However,  because  the  time  regimes  of  interest, 
though  less  than  the  closing  times,  are  slow  compared  to  electron  collision 
frequencies,  an  equilibrium  or  quasi-equilibrium  kinetic  treatment  is  possible. 

We  have  recently  built  a  kinetic  model  of  air  breakdown  capable  of  predicting 
spark  gap  closing  times  Tq  for  a  both  over—  and  under— volted  spark  gaps  [Ref.  2]. 
The  closing  time  is  determined  primarily  by  ohmic  heating  which  eventually 
enables  thermal  ionization  to  take  over.  EMP  pulse— widths  may  be  long  enough 
to  permit  sufficient  ohmic  heating  to  cause  breakdown.  We  have  found  that 
under— volted  gap  closure  is  explained  by  electron  detachment  from  0*.  Simulated 
closing  times  agree  with  empirical  data  in  a  range  from  10  »  over— volted  to  10  * 
under— volted,  which  varies  approximately  as; 

pTc  =  7.4»10H  (E/p)-3-4< 

[Ref.  3]  where  p  =  1.29  kg/m^  in  1  atm  air,  and  E  is  in  V/m. 

We  set  out  to  build  a  quasi-equilibrium  kinetic  model  of  HF  air  breakdown  based 
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on  onr  spark  breakdown  model.  An  understanding  of  the  time— dependent 
breakdown  process  is  the  most  critical  prerequisite  to  building  a  complete 
phenomenological  model  of  HF  arcing  in  apertures  illuminated  by  BMP.  The 
BMP  energy  leaking  through  an  aperture  will  depend  primarily  on  how  long  it 

takes  for  the  arc  to  form  ("100  ns);  and  to  a  lesser  degree  on  how  well  the  arc 

plasma  shields  against  further  energy  transfer  after  it  is  well— formed.  Furthermore, 
the  fast  collapse  of  the  voltage  across  the  aperture  can  create  very  high  frequency 
noise  inside  the  shielded  enclosure,  thus  actually  enhancing  the  high-frequency 

content  of  the  transferred  BMP.  The  Phase  I  program  reduced  this  highest  risk 
factor  by  examining  the  critical  phenomenology  first,  and  demonstrating  the 
adequacy  of  an  innovative  modeling  methodology  to  model  HF  breakdown  processes. 

B.  PHASB  I  OBJBCTIVBS. 

The  ultimate  goal  is  to  develop  and  verify  a  self— consistent  theoretical  model  of 

air  breakdown  in  apertures  illuminated  by  BMP.  That  goal  was  supported  by  the 
following  Phase  I  objectives; 

1.  Determine  streamer  behavior  in  feasible  B— pulse  profiles.  Key  output 

parameters  are  the  speed  of  propagation  and  nominal  electron  density  created  by  a 
streamer. 

2.  Determine  the  adequacy  of  a  quasi— equilibrium  kinetic  model.  Key 

parameters  are  the  characteristic  times  for  reaching  equilibrium  values  of  ionization, 
attachment  and  other  kinetic  rates,  as  they  compare  to  pulse— width. 

3.  Plan  for  complete  model  development.  Phase  I  0— D  and  1— D  models  would 

uncover  physics  and  numerical  issues  that  will  be  more  challenging  in  higher 
dimensions.  We  would  determine  the  appropriate  dimensionality  (1— D,  2— D  or 

3— D),  and  develop  a  plan  for  developing  such  a  model  in  Phase  II. 

4.  Plan  for  experimental  verification.  Phase  I  analysis  would  identify 
uncertainties  in  the  model.  We  would  develop  detailed  plans  to  remove  or  reduce 
the  most  significant  uncertainties  with  experiments  which  can  be  performed  within 
the  scope  of  Phase  II. 
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Those  objectives  were  supported  by  five  technical  tasks: 


1.  Boltzmann  Analysis.  We  used  a  fully  time— dependent  numerical  solution  of 
the  Boltzmann  operation  to  quantify  the  relaxation  times  of  the  electron  transport 
coefficients. 

2.  Streamer  Modeling.  We  implemented  the  Quasi— Equilibrium  Methodology 
(QEM)  in  a  single  1— dimensional  electric  field  solving  code  (elfid)  to  determine 
the  effect  of  QEM  on  streamer  propagation. 

3.  Kinetic  Modeling.  We  analyzed  the  significant  physics  dominating 

breakdown  times  in  a  0— dimensional  kinetics  code  to  determine  feasible 
simplifications. 

4.  Advanced  Model  Plan.  We  identified  the  dominant  phenomenology,  and 
included  a  complete  plan  for  developing  a  self-consistent  breakdown  model  in  a 
Phase  II  Program  Plan. 

5.  Experimental  Verification  Plan.  We  included  a  detailed  plan  for  a  Phase  II 
experiment  to  verify  the  Phase  II  model  in  the  Phase  II  Program  Plan. 

C.  OVERVIEW. 

Two  main  issues  were  examined  in  Phase  I.  The  first  issue  is  non— equilibrium 
transport;  the  second  is  the  breakdown  kinetics. 

A  clear  phenomenological  description  of  the  electrical  breakdown  of  high-pressure 
gaseous  insulators  has  emerged  only  in  the  last  decade  or  so.  Typically  a 
micro— irregularity  (accidental  or  intentionally  triggered)  will  create  a  sufficient  local 
disturbance  to  launch  a  streamer.  By  "streamer"  we  mean  a  glowing  ionization 
wave  which  enhances  the  E— field  ahead  of  its  tip,  creating  intense  local  excitation 
and  ionization  (about  10 e/cm^).  Streamers  are  usually  filamentary,  because  the 
wave  velocity  (up  to  lO^  cm/s)  increases  with  decreasing  radius  of  curvature, 
making  the  tip  self-sharpening.  Once  the  filamentary  streamer  connects  to  the 
other  electrode,  the  partiaUy  ionized  filament  draws  current;  not  enough  to  cause 
voltage  collapse  for  reasonably  large  circuit  capacitance,  but  enough  to  cause 
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significant  ohmic  heating.  When  the  filament  temperature  exceeds  about  10,000  K, 
thermal  ionization  dominates.  At  that  point,  the  system  becomes  strongly 
non-linear,  because  the  higher  the  temperature,  the  faster  the  thermal  ionization; 
the  higher  the  ionization,  the  faster  the  ohmic  heating.  The  rapidly  rising  current 
causes  sudden  voltage  collapse,  and  the  filamentary  glow  discharge  then  quickly 
evolves  into  an  arc. 

We  recognized  that  the  fast  streamer  time  scales  (significant  changes  in  local 
E— field  and  electron  density  in  a  few  picoseconds)  were  comparable  to  collision 
frequencies  (a  few  per  ps).  Therefore,  equilibrium  fluid  transport  techniques  are 
questionable.  However,  since  the  deviation  from  equilibrium  are  not  likely  to  be 
very  great,  a  perturbation  method  may  adequately  model  this  quasi— equilibrium 
regime.  This  would  avoid  expensive  particle— pushing  techniques  (e.g. 
Monte— Carlo),  which  are  more  appropriate  when  the  transport  is  nearly  ballistic 
(collissionless)  than  when  it’s  collision— dominated.  Thus,  the  objective  of  our 
Quasi— Equilibrium  Methodology  was  to  extend  fluid  model  techniques  with 
first— order  non— equilibrium  effects. 

Filamentary  streamers  may  not  form  in  the  small  micro— gaps  in  the  aperatures  of 
EMP-hardened  structures,  because  their  sizes  are  comparable  to  expected  streamer 
diameters.  However,  the  electrons  emitted  from  a  cathode  are  never  in  equilibrium 
with  the  local  electric  field  in  the  gas,  and  the  disturbance  due  to  the 
discontinuities  in  E— field,  and  the  density,  initial  velocity  and  mean  energy  of  the 
injected  electrons  take  a  certain  distance  to  relax  to  equilibrium  in  the  gas 
interior.  This  distrubance,  which  persists  in  steady-state,  is  known  as  a  cathode 
fall.  Viewed  from  the  perspective  of  a  fluid  "particle"  (a  clump  of  electrons),  the 
cathode  boundary  conditions  are  exactly  analogous  to  temporal  initial  conditions. 
Whether  the  cathode  fall  in  a  micro— gap  with  increasing  voltage  evolves  into  a 
Townsend  avalanche  or  an  ionization  wave  (planar  or  filamentary),  it  is  clear  that 
non— equilibrium  effects  must  be  considered.  Our  Quasi— Equilibrium  Methodology 
(QEM)  emphasizes  the  temporal  approach. 

Our  methodology  (QEM)  is  based  on  assuming  that  the  transport  coefficients  relax 
toward  their  equilibrium  values  on  quantifiable  time  scales.  The  implications  of 
non-equilibrium  effects  are  examined  in  detail  based  on  or  der-of— magnitude 
estimates  in  Chapter  II.  One  interesting  consequence  is  an  electron  velocity  over— 
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shoot  for  sudden  El— field  changes,  which  can  significantly  affect  streamer 
propagation.  The  quantification  of  relaxation  times  based  on  time— dependent 
Boltzmann  equation  solutions  is  described  in  detail  in  Chapter  III.  The  main 
conclusion  is  that  the  transport  coefficients  relax  on  time  scales  comparable  to  the 
energy  transfer  time.  The  non-equilibrium  effects  are  important  in  modeling  fast 
streamers  or  planar  ionization  waves  which  trigger  the  breakdown  process,  or  in 
modeling  the  cathode  fall  or  Townsend  avalanche  steady— state  which  defines  initial 
conditions  for  the  kinetic  breakdown  model.  The  effect  on  streamer  propagation  is 
shown  in  Chapter  IV. 

The  breakdown  process  itself  occurs  on  a  time  scale  long  compared  to  relaxation 
times.  Chapter  V  describes  our  assessment  of  the  potential  to  simplify  the  full 
kinetic  model  of  air  breakdown,  in  order  to  make  it  practical  for  1— dimensional  or 
multi— dimensional  codes.  The  main  conclusion  is  that  factors  of  2  to  3  in 
computational  speed  can  be  gained.  Chapter  VI  discusses  various  issues  involved 
in  developing  a  modeling  strategy.  It  develops  in  detail  the  rationale  for  our 
phenomenological  approach  for  Phase  11. 

The  modeling  approach  is  included  in  the  Phase  II  Program  Plan,  which 
constitutes  Chapter  VII.  It  includes  the  model  development  proper,  validation 
procedures,  and  implementation  in  an  existing  EMP  coupling  code.  The  validation 
includes  the  verification  experiment  plan.  The  demonstration  EMP  code  is  the 
POLYANA  code  developed  by  BDM,  and  will  require  their  participation  in 
Phase  II.  Phase  III  will  include  marketing  of  the  QEM  module  as  a  stand-alone 
product. 

Appendices  cover  two  esoteric  topics.  Appendix  A  describes  the  derivation  of  the 
Relaxation  to  a  Sliding  Quasi-Equilibrium  (RSQE)  algorithm  which  makes  the 
QEM  approach  computationally  affordable.  Appendix  B  considers  the  effect  of 
dust  particles  on  the  kinetics  and  transport,  and  demonstrates  it  can  be  neglected. 
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n.  NON-EQUILIBRIUM  ELECTRON  TRANSPORT. 


The  key  to  the  Quasi— Equilibrium  Methodology  (QEM)  is  the  time-dependent 
Boltzmann  equation  solution.  However,  some  of  the  expected  behavior  can  be 
examined  with  order-of— magnitude  estimates  based  on  reasonable  assumptions  and 
approximations.  This  chapter  examines  the  expected  non-equilibrium  electron 
transport  behavior,  specifically  electron  velocity  and  electron  multiplication. 

Since  the  momentum  transfer  process  is  dominated  by  the  lower  energies  in  the 
electron  energy  distribution  fimction  (EEDF),  we  should  expect  the  relaxation  of 
the  momentum  transfer  collision  frequency  to  occur  on  the  same  time  scale  as 
energy  transfer.  An  estimate  of  the  energy  transfer  time  can  be  obtained  as 
follows.  At  equilibrium,  the  energy  losses  must  equal  the  gains  from  ohmic 
heating.  Therefore,  the  energy  transfer  time  is  the  ratio  of  mean  electron  energy 
<  S>  to  the  ohmic  heating  rate.  For  3  MV/m  in  1  atm  air,  the  mean  electron 
energy  is  about  3  eV.  Thus,  the  relaxation  time  is  about 


r 


<S> 
e  E  Vd 


10-iis 


where  e  is  the  elementary  charge.  In  comparison,  the  equilibrium  drift  velocity, 
Vd,  for  a  nominal  E  =  3  MV/m  in  1  atm  air  is  about  10®  m/s.  Thus,  the 
momentum  transfer  collision  frequency  is  about 


1/  =  ®  s  5*1012/8 
m  Vd  ' 

where  m  is  the  electron  mass.  Therefore,  we  should  expect  the  momentum 
transfer  relaxation  time  r  to  be  about  50  collision  times  (l/i^). 

To  illustrate  what  this  means,  consider  a  homogeneous  E— field  turned  on 
"instantaneously"  Jast  compared  to  the  collision  frequency).  Neglecting  spatial 
terms,  this  is  an  initial— value  problem,  requiring  integration  only  in  time.  Before 
the  step,  all  values  are  in  equilibrium  for  E  =  0.  The  initial  velocity  is  0,  but 
the  initial  collision  frequency  is  not  0,  because  random  motion  produces  collisions 
even  without  an  L— field.  From  low  air  mobility  data  [Ref.  4]  we  can  assume  an 
initial  value  i/(0)  a  1/ps  (vs  the  equilibrium  value  i/q  =  5/ps).  With  those 
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conditions,  we  solved  the  Ordinary  Differential  Equation  (ODE)  pair 

V  =  ^E-V:-  aad  i-  =  (Kq-i/j/r 

for  time-dependent  (non-equilibrium)  electron  velocity,  V,  and  collision  frequency, 
u,  by  the  Relaxation  to  a  Sliding  Quasi— Equilibrium  (RSQE)  method 
[Appendix  A].  Figure  1  shows  the  collision  frequency  relaxing  toward  the 

equilibrium  value  of  5/ps.  Notice  there  is  still  a  significant  gap  after  10  ps, 

which  is  one  relaxation  time  r.  Figure  2  shows  the  behavior  of  the  electron 

velocity.  Notice  V  relaxes  to  the  Sliding  Quasi— Equilibrium  (SQE)  curve  much 

faster  than  the  SQE  relaxes  to  the  eventual  equilibrium  drift  velocity  Vd-  The 
test  case  conditions  instantaneously  increase  the  driving  term  E,  but  the 
drag— inducing  collision  frequency,  u,  increases  much  more  slowly.  The  result  is 
that  the  electron  fluid  velocity  overshoots  the  equilibrium  value  Vd  and  then 
relaxes  toward  it  from  the  high  side  on  the  energy  transfer  time  scale. 

This  phenomenon  can  have  a  profound  effect  on  the  behavior  of  the  streamer 
propagation  model  and  a  non-equilibrium  cathode  fall  model.  The  tendency  of  the 
electron  fluid  to  "over-react"  to  the  sudden  change  in  E— field  will  make  the 
streamer  tip  sharper  (more  shock— like),  and  may  explain  why  streamer  models 
based  on  equilibrium  transport  generally  underestimate  the  streamer  propagation 
velocity.  Consider  that  the  integral  of  the  electron  velocity  is  a  total 
displacement,  and  the  charge  separation  is  what  distorts  the  field  at  an  ionization 
wave— front.  The  total  displacement  J  V  dt  computed  by  the  RSQE  method  out 
to  several  r— values  is  about  1.4*  10‘®  m  more  than  using  the  equilibrium  Vd.  For 
a  nominal  lO^o  e/m^  (lO^^/cm^),  such  an  error  in  charge  separation  will  translate 
into  an  error  in  computed  space— charge  field  of 

e  Ne/fo  «  2.6  MV/m 

where  Ng  is  electron  density,  and  Cq  is  the  permittivity  of  free  space.  That  is 
not  a  negligible  effect,  even  though  the  "overshoot"  occurs  on  a  time— scale  of 
about  10  ps. 

A  reasonable  upper  bound  for  both  an  ionization  wave-front  and  for  a  typical 
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cathode  fall  (CF)  in  1  atm  air  is  about  10  MV/m.  For  that  field,  the  drift 
velocity  increases  to  about  2*105  m/s,  leaving  the  collision  frequency  about  the 
same  as  assumed  above  for  3  MV/m.  Furthermore,  the  mean  energy  increases  to 
about  5  eV,  leaving  the  energy  transfer  time  about  the  same.  Therefore,  the 
relaxation  time  scales  of  the  above  analysis  are  insensitive  to  E— field,  although  the 
equilibrium  drift  velocity  Vd  increases  approximately  in  proportion  to  E.  The 

total  displacement  error  SX  if  non-equilibrium  effects  are  ignored  estimated  above 

for  3  MV/m  would  double  (proportional  to  Vd)  for  10  MV/m,  and  so  would  the 
resulting  space— charge  field  error  estimate  6E. 

The  effect  of  non— equilibrium  ionization  is  harder  to  evaluate.  The  field 

distribution  in  both  an  ionization  wave  and  a  cathode  fall  is  dominated  by  the 
space-charge  effect,  and  that  is  dominated  by  the  electron  multiplication  (given  by 
the  ionization  rate  coefficient).  The  ionization  rate  is  dominated  by  the  high- 

energy  tail  of  the  EEDF,  which  is  expected  to  relax  on  a  slower  time  scale  than 

the  mean  energy.  For  10  MV/m  in  1  atm  air,  the  equilibrium  ionization  rate  Ui 
is  about  3)‘10Vs-  This  leads  to  CF  scale— size  estimates  in  the  regime  of  10*<  m 
(10 '2  cm),  which  is  comparable  to  the  expected  micro— crack  sizes.  However,  the 
electron  emission  mechanisms  are  unlikely  to  produce  electrons  above  the  ionization 
thresholds  (>12  eV),  but  will  produce  plenty  above  the  dissociative  attachment 
threshold  ('3  eV).  Therefore,  the  appropriate  initial  (boundary  condition) 
ionization  rate  J/i(0)  is  much  smaller  than  the  attachment  rate  i/a(0),  which  will 
be  close  to  its  equilibrium  value  of  about  10®/s.  If  the  ionization  rate  relaxation 
time  is  much  slower  than  l/i/j,  the  electron  multiplication  in  a  sub— mm  aperture 
will  be  much  less  than  calculated  by  equilibrium  models. 

For  illustrative  purposes,  we  solved  the  ODE  pair 

I^e  =  s  —  Ne(t'a-i^i)  and  i'i  =  (t'q-i'i)/ri 

with  the  initial  conditions  Ne(0)  =  lO^/cm^  and  =  0,  a  constant  equilibrium 

ionization  rate  Uq  =  3/ns,  and  relaxation  time  ri  varying  from  0.1  to  3.0  ns. 
The  background  source  S  was  chosen  so  the  pre-existing  steady-state  value  S/i/a 
matched  the  initial  electron  density  Ne(0).  The  ionization  rate  grows  quickly,  so 
the  net  (i^a— t^i)  becomes  a  gain  (multiplication)  term  almost  immediately. 
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Figure  3  shows  the  relaxation  of  the  net  gain  rate  (i^i— t'a),  which  is  positive 
{ui  >  I/a)  most  of  the  time.  Figure  4  shows  the  resulting  electron  multiplication 
curves.  Note  that  on  a  log-linear  plot,  positive  slope  corresponds  to  net  gain 
multiplication  rate  {ui—i/g).  The  ri  =  0.1  ns  curve  is  almost  an  equilibrium 
solution,  which  assumes  the  ionization  rate  i/j  comes  into  equilibrium  with  the 

E— field  virtually  instantaneously.  The  value  of  the  relaxation  time  n  clearly  can 
affect  the  order  of  magnitude  of  the  electron  number  density  in  a  hairline  crack 
created  by  an  overvoltage,  and  must  be  evaluated  for  the  QEM  approach. 

This  test  case  is  analogous  to  a  cathode  fall,  but  only  roughly.  In  a  real  cathode 
fall  [Ref.  5,  6,  and  7],  the  E— field  drops  off  approximately  bnearly,  but  a  nominal 
10  MV/m  is  the  right  order  of  magnitude.  The  analogous  test  case  we 

constructed  corresponds  to  a  constant  10  MV/m,  and  since  the  drift  velocity 
relaxes  to  the  equifibrium  2>«105  m/s  on  a  10— ps  time  scale,  5  ns  in  time  roughly 
corresponds  to  1  mm  in  distance. 

However,  the  E— field  is  not  going  to  be  constant,  either  in  a  CF  or  a  micro- 

crack.  Assuming  sharp  edges  enhance  the  peak  field,  since  the  equiUbrium 
ionization  rate  is  very  sensitive  to  E— field,  it  is  not  unrealistic  to  assume  the 
equilibrium  rate  ui  starts  at  about  10/ns,  but  decreases  to  much  less  than  the 

nominal  3/ns.  Therefore,  we  repeated  the  CF— analogue  analysis  with  the 

equilibrium  ui  decreasing  from  1C  to  0.3  /ns.  Figure  5  shows  the  resulting 

instantaneous  ionization  rates  curves.  The  increase  for  slower  ti  at  late  times 
reflects  the  fact  that  slower  relaxation  times  delay  the  cooling  of  the  high-energy 
tail  of  the  EEDF  as  much  as  its  heating  at  early  times.  Figure  6  shows  the 

resulting  electron  multiplication  curves.  Notice  the  lower— rj  curves  eventually 
exceed  lO'^/cm^.  Near  those  levels,  streamer— like  ionization  waves  may  form  and 
ionize  the  whole  gap. 

For  0.1— mm  micro— gaps,  the  avalanche  discharge  would  be  truncated  at  a  position 
corresponding  to  about  0.5  ns.  However,  the  effective  heating  and  electron 

multiplication  would  accumulate  for  a  sinusoidal  driving  signal  from  one  half— cycle 

to  the  next.  To  first  order,  each  0.5— ns  interval  represents  a  half— cycle. 

Considering  that  the  polarity  changes  every  half-cycle,  the  new  cathode  initial 
conditions  are  the  previous  half-cycle’s  final  conditions  at  the  anode.  If  we 
assume  breakdown  occurs  in  the  cycle  where  Ne  reaches  10i</cm3,  we  conclude 
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that  several  cycles  will  be  needed. 


Non— equilibrium  effects  on  the  electron  multiplication  will  be  negligible  if  the 
ionization  coefficient  relaxation  time  ri  is  slow  compared  to  the  multiplication  time 
1/i/i,  but  will  dominate  if  ri  >  Ijui-  The  time— dependent  Boltzmann  analysis 
leads  to  quantitative  estimates  of  all  relaxation  times. 
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Figure  2.  Relaxation  of  V:  0  Vd(105m/s)  as  1/v. 
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Figure  3.  Net  Gain— Loss  Rates  for  Constant-E  Case. 
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Figure  5.  Net  Gain— Loss  Rates  for  Deaeasing— £  Case. 
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m.  TIME-DEPENDENT  BOLTZMANN  ANALYSIS  OF  KINETIC 
RELAXATION. 

The  basis  of  the  Quasi— Equilibrium  Methodology  being  developed  by  Tetra  is  a 
quantitative  characterization  of  the  relaxation  of  the  electron  transport  and  kinetic 
coefficients  for  fast  changes  in  local  electric  field,  E.  Such  fast-changing  E— fields 
are  expected  in  avalanche— driven  ionization  waves,  whether  filamentary  streamers  or 
planar  micro— gap  breakdown  discharges.  This  chapter  documents  the  calculation  of 
characteristic  relaxation  times  for  1— atm  dry  air  from  solutions  of  the  time- 
dependent  Boltzmann  equation. 

The  integro— differential  Boltzmann  equation  is: 

If  *  r,  .  5  E.?v]  f(r,v,t)  = 

The  operators  in  the  first  parentheses  are  the  time  derivative,  the  gradient  in  real 
space,  and  the  inner  product  of  the  acceleration  vector  with  the  gradient  in 
velocity  space,  respectively.  The  operand,  f,  is  the  electron  distribution  function, 
and  the  right-hand-side  is  the  total  effect  of  collisions.  The  balance  between  all 
heating  and  cooling  mechanisms  establishes  the  equilibrium  distribution  function. 
In  the  absence  of  quantized  energy— selective  processes,  the  steady— state  solution  is 
a  Maxwellian  distribution,  proportional  to  ea:p(— ^/kTg).  However,  energy-selective 
collisions  distort  the  distribution  function  from  a  Maxwellian. 

If  the  E— field  varies  on  a  larger  scale— size  than  the  distribution  function,  then  the 
Local  Field  Approximation  (LFA)  applies,  and  the  Vr  operator  may  be  neglected. 
We  shall  assume  LFA  applies  either  directly  or  through  the  use  of  the  method  of 
characteristics,  which  by  following  electron  fluid  particles  (at  the  mean  directed 
velocity),  turns  spatial  variations  into  temporal  variations.  Any  coefficient  Kp  may 
be  calculated  by  convoluting  the  instantaneous  f  with  the  cross  section  for  that 

process  for  which  the  cross  section  is  Cp  as  follows: 

V  /  (Tpf  </2  ^/m  dg 

=  /  f 

If  the  time  scale  on  which  the  E— field  varies  is  slow  enough,  then  the  steady— 
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state  EEDF  may  be  used  to  calculate  equilibrium  transport  and  kinetic  rate 
coefficients. 

The  moments  of  the  Boltzmann  equation  take  the  form  of  energy  and  momentum 
conservation  and  particle  continuity  equations.  These  can  be  solved  using  the 
equilibrium  transport  and  kinetic  rate  coefficients,  together  with  the  coupled  E— field 
equations  to  form  a  fluid  model  of  a  gas  discharge.  The  Quasi— Equilibrium 
Methodology  extends  the  fluid  approach  into  the  non— equilibrium  domain  by 
approximating  the  relaxation  of  the  transport  and  kinetic  coefficients  by  use  of 
characteristic  relaxation  times.  The  purpose  of  our  time— dependent  Boltzmann 
analysis  is  to  quantify  the  relaxation  times  for  limited  deviations  from  the 
equilibrium  conditions. 

For  the  analysis,  we  use  the  time-dependent  Boltzmann  equation  solving  code 

ELENDIF  [Ref.  8],  which  uses  implicit  time  integration  techniques  to  insure  stable 
relaxation  to  the  steady-state  solution.  The  electron— impact  cross  section  data  are 
from  [Ref.  9]  for  Nj,  [Ref.  10]  for  O2,  [Ref.  11]  for  Ar,  and  [Ref.  12]  for  HjO. 
The  time— dependent  relaxation  calculation  includes  superelastic  collisions, 
attachment,  recombination  and  ionization.  We  modified  the  code  to  use  a 

geometric  time  progression,  and  to  calculate  an  output  at  each  time  step:  the 

non— equilibrium  electron  energy  distribution  function  (EEDF),  the  transport 
coefficients,  the  energy  flow  rates,  and  the  inelastic  and  superelastic  rate 
coefficients. 

The  relaxation  times  can  best  be  calculated  for  a  step— function,  instantaneous 
change  in  El— field,  with  initial  conditions  not  too  far  from  the  new  equilibrium. 
Therefore,  we  stepped  up  in  E/N  in  increments  of  factors  of  2  or  so.  The  EEDF 
was  initialized  to  a  0.1-eV  Maxwellian,  and  the  E/N  was  first  set  to  20  Td 
(20*10'i^  V*cm2).  Convergence  was  achieved  in  less  than  1  ns,  but  the  time- 

sequence  continued  for  several  ns.  The  E/N  was  then  stepped  up  to  50  Td,  using 
the  20— Td  steady— state  solution  as  the  initial  EEDF  for  the  second  geometric  time 
sequence  (0.1  ps  to  4.6  ns).  This  continued  up  the  E/N  sequence  20,  50,  100, 
200,  500  Td,  and  then  we  stepped  down  200,  100,  50  and  20  Td.  The  resulting 
time-dependent  EEDF  sequence  for  the  up— stepping  half  is  shown  in  Figure  7. 

The  structure  between  6  and  8  eV  at  very  early  times  is  due  to  superelastic 
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collisions  with  a  nominal  population  of  10-<  N2(A)  (i.e.,  one  molecule  in  10,000  is 
assumed  in  the  exdted  state).  Other  runs  with  10'^  and  lO'^  N2(A)  showed 
the  peak  is  indeed  sensitive  to  the  assumed  excitation  level,  affecting  the  rate 
coefficients,  especially  for  processes  with  thresholds  in  that  vicinity,  such  as  O2 
ionization.  However,  the  relaxation  behavior  (discussed  in  detail  below)  is  virtually 
unaffected.  Nonetheless,  it  would  be  wise  for  a  complete  breakdown  model  to 
account  for  the  dependence  of  rate  coefficients  on  the  N2(A)  population. 

We  expected  the  high-energy  tail  to  relax  to  a  steady  state  on  a  slower  time 
scale  than  the  low-energy  portion.  It  seems  the  normalization  of  the  EEDF 
causes  changes  at  all  energies  at  all  times.  This  agrees  qualitatively  with  E.  E. 
Kunhardt’s  results  [Ref.  13]  in  pure  N2.  He  solved  the  time— dependent  Boltzmann 
equation  with  a  totally  different  code,  and  observed  that  the  EEDFs  for  different 
initial  conditions  relaxed  to  a  common  EEDF  on  the  energy  relaxation  time  scale, 
and  then  to  the  final  steady— state  on  a  slower  scale. 

A  common  approach  to  modeling  non-equilibrium  transport  is  that  used  by 
P.  Bayle,  J.  Vaquie  and  M.  Bayle  [Ref.  14],  which  we  shall  call  the  BVB  model. 
They  assume  all  transport  and  rate  coefficients  are  functions  of  electron 
temperature  or  mean  energy  <^>.  By  solving  the  energy  conservation  equation, 
they  calculate  the  non— equilibrium  mean  energy,  and  estimate  the  rate  coefficients 
by  interpolation  as  functions  of  that  mean  energy.  That  is  reasonable  if  each 
non-equilibrium  EEDF  approximates  an  equilibrium  shape  at  an  intermediate  E/N 
value.  If  that  hypothesis  is  true,  then  the  non-equilibrium  coefficients  we 
calculate  by  the  time— dependent  Boltzmann  analysis  should  be  smooth  functions  of 
<S>.  That  suggests  a  simple  graphical  test  of  the  hypothesis. 

Figure  8  shows  the  time-dependent  behavior  of  mean  energy  in  our  Boltzmann 
simulation.  Since  it  is  monotonic  in  each  half  (up  and  down  in  E/N),  we  can 
substitute  <  ^>  for  time  as  the  independent  parameter.  Figure  9  shows  reduced 

collision  frequencies  ://N  vs  <^>.  Notice  that  energy  transfer  is  consistently  over 
one  order  of  magnitude  slower  than  momentum  transfer.  It  stands  to  reason  that 
energy  transfer,  which  depends  on  inelastic  collisions,  should  be  slower  than 

momentum  transfer,  which  is  dominated  by  the  more  frequent  elastic  collisions. 
The  lack  of  large  deviations  is  encouraging.  It  explains  why  the  BVB  model 

gives  reasonable  behavior;  the  difference  between  a  non— equilibrium  collision 
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frequency  and  an  interpolation  between  equilibrium  values  is  never  large.  However, 
the  systematic  trend  of  the  error,  more  obvious  in  energy  transfer  than  in 
momentum  transfer,  is  disturbing.  The  "orbits"  reminiscent  of  hysteresis  loops 
appear  consistently  in  plots  of  all  rate  coefficients  vs  <  S>,  even  the  high  energy- 
threshold  rates  shown  in  Figure  10.  Close  examination  shows  the  rate  coefficients 
are  clearly  relaxing  faster  than  <  S>.  Thus,  the  slow  relaxation  of  <  S>  (not  to 
be  confused  with  the  relaxation  of  the  energy  transfer  rate)  casts  doubt  on  the 
basic  premise  of  the  BVB  approach. 

A  closer  examination  of  details  of  the  Boltzmann  output  (specifically  the  energy 
balance)  revealed  an  explanation.  The  instantaneous  net  energy  rate  can  be 

expressed  as  follows: 


<#>  =  -  E  Ve  -  <g> 

m  ^  e 

The  heating  term  E  Vg  e/m  is  changing  on  the  same  time  scale  as  the  energy 
transfer  collisional  frequency  i/ g  Consequently,  the  two  nearly  balance,  so  the 
mean  energy  <^>  changes  much  more  slowly  than 

We  propose  an  alternate  methodology  based  on  quantifying  a  characteristic 
relaxation  time.  Figure  11  shows  how  the  time— dependent  elastic  collision 
frequencies  relax  toward  equilibrium.  Most  of  them  show  close  to  an  exponential 
decay  of  the  difference  |K— Kq|  (a  straight  line  on  a  log-linear  plot).  The 
desired  relaxation  times,  r,  must  be  extracted  from  the  temporal  behavior  of  the 
calculated  coefficients,  K.  Assuming  the  governing  equation  is  of  the  form 


K  =  (K,-K)/r 

where  Kq  is  the  equilibrium  value,  the  average  relaxation  time  in  a  finite  time 
step  At  can  be  estimated  as 


r 


At  / 


In 


fKqzdLl 

iKq-K’J 


The  results  of  that  calculation  for  momentum  transfer,  which  we  expected  to  relax 
on  a  time  scale  comparable  to  energy  transfer,  and  for  ionization,  which  we 
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expected  to  relax  on  much  slower  time  scales,  is  shown  in  Figure  12. 

The  algorithm  used  is  cleaxly  inaccurate  at  late  times  when  |Kq— K|  <<  Kq,  and 
differencing  errors  make  the  results  meaningless.  Calculations  at  those  times  are 
not  included  in  Figure  12.  Nonetheless,  the  behavior  of  the  calculated  r  values  is 
rather  peculiar.  The  calculations  at  very  early  times  are  not  particularly 
meaningful,  either.  The  algorithm  is  most  accurate  when  At  «  r.  The  "r  =  t 
slant  lines"  in  Figure  12  are  near  that  optimal  point.  If  we  use  the  intersecting 
points  (within  the  box  for  each  E/N  value)  as  "nominal"  r  values,  we  find  that 
the  high-energy  processes  are  barely  slower  to  relax  than  momentum  transfer,  and 
both  are  about  10  ps,  which  is  just  about  what  we  expected  for  the  energy 
transfer  relaxation  time,  as  shown  in  Chapter  II. 

Notice  that  the  ionization  r— curves  consistently  approach  the  nominal  intersection 
point  from  above,  and  the  momentum  transfer  r— curves  consistently  approach  it 
from  below.  This  consistent  pattern  indicates  a  physically  significant  difference  in 
the  relaxation  behaviors.  However,  it  is  not  clear  how  to  exploit  this  observation. 

There  are  two  main  conclusions;  1)  The  momentum  transfer  collision  frequency 
relaxes  on  about  the  energy  transfer  time  scale  ('10  ps).  2)  The  high  energy- 
threshold  rates  relax  on  too  fast  a  time  scale  (corapa^ed  to  the  rates  themselves) 
for  non— equilibrium  effects  to  be  of  much  consequence.  Therefore,  the  Quasi- 
Equilibrium  Methodology  can  be  based  on  the  assumption  that  each  coefficient  K 
relaxes  to  its  equilibrium  value  Kq  as: 


K  =  (K,-K)/t 


If  a  detailed  Boltzmann  analysis  is  unavailable,  one  may  use,  to  first  order,  the 
energy  relaxation  time: 

~  1  _  <  ^ 
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Figure  8.  Non— Eq.  Mean  Energy  for  Successive  Relaxations. 
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Figure  9.  Colliaon  Frequency  Orbits  in  Energy  Space. 
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Figure  10.  Selected  Kinetic  Rates  Orbits  in  Energy  Space. 
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Figure  11.  Elastic  Collision  Relaxations  in  Time. 
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Figure  12.  Calculated  Relaxation  Times. 
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Time  =  0.1  ps  to  4.6  ns  for  each  E/N 


IV.  NON-EQUILIBRIUM  EFFECTS  ON  STREAMER  PROPAGATION. 

The  Quasi— Equilibrium  Methodology  (QEM)  purports  to  model  non— equilibrium 
effects  in  fast  ionization  waves.  A  cogent  argument  was  made  in  Chapter  11  that 
non— equilibrium  effects,  especially  the  velocity  overshoot  effect,  would  have  a 
profound  effect  on  streamer  propagation.  The  purpose  of  this  chapter  is  to 
document  a  numerical  test  of  this  hypothesis  in  the  1— dimensional  electric  field 
solver  ELFID. 

The  ELFID  code  solves  the  current  continuity  equation 

V-(J+D)  =  0 

where  J  =  eNgVe  is  conduction  current,  and  D  =  e  E  is  the  displacement  vector, 
whose  time  derivative  is  displacement  current.  The  ELFID  code  includes  a  gas 
kinetics  module  which  updates  the  conductor  density  eNg  and  velocity  Vg  over  a 
finite  time— step  At,  given  initial  and  final  E— field  arrays. 

The  code  was  made  operational  on  an  80386-based  PC  using  the  Lahey  Fortran 
compiler  to  produce  protected-mode  code,  thus  avoiding  the  memory  limitations  of 
DOS.  The  gas  kinetics  module  was  modified  to  use  the  QEM  current  integration 
algorithm,  coded  in  subroutine  CURINT.  That  routine  integrates  the  two  pairs  of 
ordinary  differential  equations  (ODE) 

Ng  =  S  -  R  Ng  R  =  (Vaq-R)/r 

V,  =  ±  E  -  r  Ve 

where  S  is  an  external  ionization  source;  the  loss  rate,  R,  tends  to  the  net 
difference  between  equilibrium  attachment  (/?q),  and  Townsend  ionization  (oq)  rates; 
the  collision  frequency,  u,  tends  toward  the  equilibrium  value,  i/q.  The  relaxation 
times,  r,  for  both  ODE  pairs  were  set  to  the  energy  transfer  time: 

< 

’■  =  STTVT 
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In  order  to  determine  the  effect  of  the  QEM  corrections,  it  was  necessary  to 
artificially  defeat  them.  For  1  atm  air  near  3  MV/m,  the  energy  transfer  time  is 
about  lO'ii  s,  but  the  momentum  transfer  frequency  is  about  5*10i2/s.  Therefore, 
by  artificially  reducing  the  relaxation  time  by  a  factor  of  10  *3,  we  made  the 
relaxation  of  transport  coefficients  virtually  instantaneous. 

The  test  case  chosen  was  a  spherical  ionization  wave  driven  by  a  rising  voltage. 
The  spherical  geometry  simulates  the  region  ahead  of  a  filamentary  streamer  tip, 
and  insures  a  field  enhancement  sufficient  to  launch  an  ionization  wave.  A  similar 
simulation  in  SFs  is  reported  by  R.  Morrow  [Ref.  15]  with  the  same  essential 
qualitative  features.  While  both  cases  are  better  described  as  voltage— driven 
growing  coronas,  they  contain  the  essential  phenomenological  features  of  freely 
propagating  streamers.  In  our  case,  the  gap  goes  from  1  mm  to  1  cm  radii, 
insuring  a  vacuum  field  enhancement  factor  of  10;  the  driving  voltage  rises  to 
30  kV  in  10  ns.  External  ionization  is  negligible,  but  the  electron  density  is 
initialized  to  IQis/m^  to  simulate  the  photoionization  ahead  of  the  glowing  streamer 
tip  [Ref.  16). 

The  results  are  compared  in  Figures  13  and  14,  each  of  which  contains  two 
"frames"  from  the  time-dependent  simulations  with  and  without  the  QEM 
correction.  The  qualitative  behaviors  are  very  similar:  the  streamers  are  launched 
at  about  the  same  time,  and  propagate  at  about  the  same  velocity.  The  electron 

avalanching  starts  when  E  exceeds  3  MV/m,  and  the  exponentially  growing 

electron  density  eventually  shunts  the  electric  field.  After  a  few  transient 

oscillations,  the  field  will  equilibrate  at  the  critical  value  of  3  MV/m.  In  effect, 
the  discrete  grid  cells  are  acting  like  glow  discharge  voltage  regulators  in  series, 
attempting  to  clamp  the  E— field  at  3  MV/m.  In  the  limit  of  very  slow  voltage 
rise  rate,  the  solution  would  be  a  series  of  steady-states  with  the  E— field  clamped 
to  3  MV/m  in  the  corona,  and  falling  off  as  1/R2  outside  the  corona.  That 
clamping  requirement  and  the  voltage  rise  rate  determine  the  "streamer"  velocity 
to  first  order. 

The  striking  difference  between  the  two  simulations  is  in  the  values  of  electron 
density  created  by  the  ionization  waves:  about  lO^i/m^  without  QEM,  but  only 

10i9/m3  with  QEM.  An  explanation  can  be  found  in  the  electron  velocity  over¬ 
shoot.  Because  the  collision  frequency  takes  many  collisions  to  relax  to 


29 


equilibrium,  the  instantaneous  velocity  overshoots  the  equilibrium  drift  velocity 

value  Vd-  This  tendency  of  the  electron  fluid  to  initially  over-react  to  the  fast- 
rising  E— field  makes  it  create  the  charge— separation  necessary  to  shunt  the  E— field 
more  quickly.  As  a  result,  fewer  electrons  are  just  as  effective  in  shunting  the 
local  over-voltage,  and  the  transition  occurs  sooner,  creating  much  lower  ionization 
levels. 

The  lower  electron  der  'ty  created  by  an  ionization  wave  can  prolong  the 
breakdown  time  significantly.  The  breakdown  time  is  roughly  the  time  it  takes 

ohmic  heating  to  raise  the  temperature  to  about  10^  K,  and  the  ohmic  heating 

rate  is  proportional  to  Ng.  Although  Ne  does  not  remain  constant  after  the 

streamer  crosses  the  gap,  the  initial  electron  density  created  by  the  streamer  does 
significantly  affect  the  total  breakdown  time.  Thus,  the  QEM  approach  to 
modeling  the  initial  phases  of  air  breakdown  is  essential. 

Another  less  obvious  difference  is  that  the  simulati^'n  without  QEM  got  bogged 

down  after  about  3.1  ns,  with  the  time  step  becoming  comparable  to  a  collision 
time  (a  few  ps).  But  the  simulation  with  QEM  continued  quite  well  beyond 

6  ns,  with  time  steps  at  worst  comparable  to  the  energy  transfer  time  (0.1  ns). 
The  later  time  behavior  is  displayed  in  Figures  15  and  16,  which  display, 
respectively,  E  and  Ng  curves  at  the  frames  nearest  the  chosen  times.  This 

ability  of  the  QEM  algorithms  to  take  reasonably  large  time  steps  will  be  very 

valuable  in  Phase  II. 
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Figiire  13.  Two  frames  from  streamer  simulation  without  QEM. 
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Figure  14.  Two  frames  from  streamer  simulation  with  QEM. 


SPHEP  TEST  with  QEM 


SPHER  TEST  with  QEM 


32 


V. 


SIMPLIFICATION  OF  THE  AIR  BREAKDOWN  KINETICS  MODEL. 


The  air  breakdown  (BD)  model  developed  by  Tetra  [Ref.  15]  tracks  the 
populations  of  16  components  with  several  dozen  reactions  among  them.  For  a 
0— dimensional  model,  computational  cost  is  not  significant,  but  for  an 
N— dimensional  code,  the  large  number  of  nodes  justifies  simplification  of  the  point— 
by— point  kinetic  computations.  This  chapter  describes  our  effort  to  simplify  the 
air  BD  kinetics  model. 


The  kinetics  of  an  electric  discharge  in  air  are  rather  complex,  and  an  over¬ 
simplified  model  is  doomed  to  a  very  narrow  domain  at  best.  Some  mechanisms 
are  obviously  necessary  in  any  air  BD  model,  including  both  electron^mpact 
phenomena  and  heavy— heavy  collisions.  We  shall  discuss  them  in  approximate 
order  of  importance,  with  the  actions  we  took  to  simplify  the  model,  including 
validation  test  runs,  leading  to  our  recommendations  for  a  simplified  model. 


It  is  clear  that  any  model  of  air  BD  must  include  the  dominant  sources  and  sinks 
of  electrons.  Those  are,  respectively,  electron  multiplication  by  molecular 
ionization: 


e' 


+  2  e* 


and  dissociative  attachment: 


e-  +  O2  -*  O-  +  0 


Electronic  excitation  followed  by  ionization  of  the  excited  species: 

e'  +  N2  e‘  +  N2(A) 
e*  +  N2(A)  N2*  +  2  e' 

was  found  to  be  important.  N2(A)  stands  for  the  A^S  state  at  6.17  eV.  The 
many  paths  by  which  electron  impacts  populate  N2(A)  can  be  combined  as  shown 
in  Figure  17  with  no  loss  of  phenomenology,  and  our  benchmark  calculations 
pr  wed  it  to  have  no  effect  on  the  model  behavior.  A  secondary  mechanism  is 
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excitation  of  the  a^E  state  at  8.40  eV  (denoted  "Nja")  followed  by  self-quenching 
ionization: 


e'  4-  N2  -♦  6'  +  N2a 
N2a  +  N2a  N2  +  N2*  +  e* 

Elimination  of  that  path  increased  the  simulated  closing  time  for  the  test  case 
(25kV/cm)  from  156  ns  to  164  ns.  The  small  gain  in  computational  efficiency  is 
not  worth  the  loss  in  accuracy. 

Because  atomic  0  and  N  react  to  produce  significant  ionization,  the  dominant 
sources  and  sinks  of  atomic  species  must  be  included.  The  significant  sources  of 
dissociation  are  impact  with  electrons,  molecules  and  atoms: 


Electronically  excited  O2  levels  above  the  dissociation  energy  immediately  dissociate. 
That  mechanism  dominates  over  direct  dissociation  of  O2,  as  shown  in  Figure  18. 
Summing  all  paths  into  one  effective  dissonation  rate  had  no  discernible  effect  on 
model  behavior.  The  atomic  population  then  produces  ionization  by  electron 

impact  or  by  association: 


1]  + 

ni 

ro2*i 

Ini  1 

^  1 

N(  ■* 

NO* 

IMj 

N2* 

The  key  to  correct  model  behavior  (why  the  gap  closes  at  all)  for  an  under- 
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volted  gap  is  molecular-impact  associative  detachment: 

+  {ni}  -  *■  +  {nIo} 

which  is  also  the  major  source  of  the  complex  molecules  O3  and  N2O  in  air 
discharges. 

The  above  processes  are  clearly  necessary,  but  they  are  not  sufficient.  Other 
processes  in  the  full  model  affect  the  BD  behavior  in  subtle  ways  not  as  easily 
understood.  Therefore,  we  sought  to  identify  negligible  processes  so  we  could 
simplify  the  model  by  elimination,  rather  than  building  up  from  what  we  know  is 
needed.  A  nominal  characteristic  life-time  for  the  "first"  component  of  each 
reaction  can  be  computed  based  on  the  maximum  population  of  the  "second" 
component  during  the  full  model  simulation.  (The  choice  of  "first"  and  "second" 
is  a  subjective  judgment.)  The  density  traces  for  the  25  kV/cm  test  case  are 
shown  in  Figure  10. 

Based  on  that  methodology,  we  determined  that  quenching  of  the  N2(A)  by  N  or 
NO  is  negligible  compared  to  self-quenching.  We  found  many  NO— impact 
reactions  could  be  eliminated  with  no  ill  effect.  The  other  eliminated  reactions 
are; 


N2  +  NO  N  +  N  +  NO  "1.^ 

N  +  NO-N  +  N  +  0  'O.lfjs 

O  +  NO^O  +  N  +  0  'O.lfjis 

0  +  NO  -*  N  +  O2  'O.bfjs 

but  we  kept  the  following  reactions: 

NO  +  (Sjl  -  N  +  0  +  |0^}  -3as 

NO  +  NO  N  +  0  +  NO  >100ns 

O2  +  NO  -»  0  +  0  +  NO  ■45ns 

NO*  +  O2  -*  NO  +  O2*  'Ins 


The  last  is  actually  a  charge— exchange  process.  Other  charge-€xchanges  were 
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neglected  because  of  their  slow  nominal  characteristic  times  even  at  10,000  K,  as 
shown: 


O  +  O2*  "*  O2  +  0  *  700^ 

O  +  NO*  ^  0*  +  NO  'lOOO/iS 

N  +  NO*  -  N*  +  NO  'SO/iS 

But  we  kept  the  following  faster  charge-exchanges: 

{n*}  +  N2  -  {§}  +  N2*  -^Qns 

and  the  "atom— exchange"  reaction 

0  +  N2  ^  N  +  NO  ■25ns 


With  the  above  eliminations  only,  the  BD  model  behavior  is  unaffected  for  the 
test  case.  Some  other  reactions  may  also  be  negligible,  but  the  gains  are 
diminishing  for  the  work  required  to  validate  the  model  after  each  elimination. 
As  it  stands,  the  simplified  air  BD  model  is  fast  enough  to  run  interactively  in 
1—2  minutes  on  a  VAX  11—780. 

The  computational  cost  is  more  sensitive  to  number  of  components  (there  is  one 
ODE  per  component)  than  to  number  of  reactions  (each  reaction  merely 
contributes  a  loss  or  gain  term  in  a  summation).  Therefore,  we  need  ways  to 
reduce  components  entirely.  The  most  promising  strategy  is  to  neglect  chemical 
composition.  That  would  reduce  the  number  of  species  (and  number  of  ODEs) 
from  16  to  the  following  9: 


M2= 


N2 

O2 

NO 


Mj 


“ISI' 


A*,  0-,  N2(A),  N2a,  M3= 


i^lo} 


The  NO  reactions  can  be  included  in  a  weighted  average  if  we  can  calculate  a 
nominal  NO  fraction.  Since  the  NO  fraction  is  significant  only  near  the  end  of 
the  test  case  simulation,  we  suggest  treating  it  as  a  function  of  temperature. 
Thus,  the  simplifiec*  model  should  reproduce  the  test  case  used  for  calibration 
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exactly.  However,  the  range  of  applicability  of  the  simplified  model  for  different 
conditions  needs  to  be  determined  with  other  test  cases.  That  complicates  the 
validation  procedure. 

The  cost  in  manpower  must  be  carefully  assessed  against  the  gains  in 
computational  speed.  The  above  method  should  produce  about  a  factor  of  2 
speed-up,  but  the  validation  procedure  requires  a  calibration  test  case  "typical"  of 
the  application,  and  suffident  exploration  of  the  multidimensional  variable  space  to 
determine  the  domain  of  applicability.  That  is  beyond  the  Phase  I  scope,  but 
may  be  fruitful  in  Phase  H. 
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Figure  18.  Dominant  Oj  Dissociation  Rates. 
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VI.  TOWARD  A  HAIRLINE  APERTURE  BREAKDOWN  MODEL. 


It  appears  that  EMP— like  radiation  may  cause  air  breakdown  in  hairline  cracks  on 
metallic  enclosures  designed  to  shield  electronic  equipment.  The  enhanced  field  in 
such  apertures  could  easily  exceed  the  DC  breakdown  level  (3  MV/m),  but  the 
decidedly  nonlinear  mechanism  is  not  well  understood.  The  purpose  of  this 
chapter  is  to  describe  an  approach  to  modeling  such  a  breakdown,  based  on  the 
Quasi— Equilibrium  Methodology  (QEM). 

It  is  difficult  to  characterize  the  apertures  accurately,  because  they  are 
unintentional.  However,  a  range  of  about  10  to  100  fjm  (lO'^— 10*^  m)  is 
reasonable.  This  is  smaller  than  the  expected  diameter  of  a  filamentary  streamer 
in  a  spark  gap  [Ref.  17,  18,  and  19].  Therefore,  we  cannot  assume  a  filamentary 
streamer  provides  the  initial  ‘lO^^*  e/cm^  which  triggers  breakdown  in  a 
conventional  (large)  air  gap  [Ref.  2]. 

The  scale  size  of  interest  is  comparable  to  the  expected  size  of  a  cathode  fall 
(CF)  at  1  atm,  and  the  electron  swarm  is  clearly  not  in  equilibrium  with  the 
local,  instantaneous  electric  field.  A  typical  CF  drops  about  100  V  across  about 
10  /xm,  for  a  nominal  field  of  about  10  MV/m.  One  of  the  best  non-equilibrium 
CF  models  reported  in  the  literature  is  that  of  P.  Bayle,  J.  Vacquie  and  M. 
Bayle  [Ref.  14],  which  we  will  refer  to  as  the  BVB  model.  The  basis  of  their 

model  is  the  assumption  that  the  non— equilibrium  transport  coefficients  can  be 
adequately  approximated  as  the  equilibrium  values  corresponding  to  the 

instantaneous  electron  temperature,  which  is  calculated  by  solving  the  energy 
conservation  equation.  They  formally  assume  a  Maxwellian  electron  energy 

distribution  function  (EEDF),  but  that  is  not  strictly  required  by  their 

methodology.  The  BVB  model  works  well  because  both  the  electron  energy  and 
transport  (momentum  transfer  collision  frequency)  are  dominated  by  the  low— energy 
portion  of  the  EEDF.  A  weakness  in  the  BVB  model  is  that  the  electron 
multiplication  (net  of  gains  by  Townsend  ionization  minus  losses  to  attachment)  is 
dominated  by  the  high-energy  tail  of  the  EEDF,  since  the  ionization  threshold  is 
12.06  eV  for  O2  [Ref.  10]  and  15.6  eV  for  N2  [Ref.  11].  Since  the  high-energy 
tail  is  populated  by  many  elastic  and  superelastic  collisions  redistributing  the 
energy,  the  relaxation  of  the  ionization  coefficient  may  occur  on  a  slower  time 
scale  than  overall  energy  transfer. 
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The  main  alternatives  we  know  of  to  handle  these  non— equilibrium  phenomena  are 
Monte  Carlo  methods  and  solution  of  the  time— dependent  Boltzmann  equation. 
Monte  Carlo  methods  are  very  robust  and  can  be  made  arbitrarily  accurate  by 
using  enough  sample  "particles",  but  they  are  computationally  very  expensive. 
Multi— temperature  models  have  been  tried  with  some  success  [Ref.  20].  These 
essentially  amount  to  solving  the  Boltzmann  equation  with  a  crude  energy— space 
grid.  The  QEM  approach  we  prefer  is  based  on  time— dependent  Boltzmann 
equation  solutions,  which  are  adequate  as  long  as  the  local  field  approximation 
holds  in  a  time— varying  sense  (or  its  equivalent  along  characteristic  lines).  At 
1  atm,  a  fluid  treatment  is  justified,  avoiding  the  expense  of  "particle— pushing". 
The  QEM  approach  we  envision  would  use  pre— computed  Boltzmann  equation 
solutions  to  characterize  the  relaxation  time  scales.  An  in-line  time— dependent 
Boltzmaim  solution  is  not  out  of  the  question,  but  it  would  be  quite  wasteful  if 
repeated  inside  the  inner  (non-linear)  loop. 

For  Phase  I,  we  did  an  order— of— magnitude  analysis  to  assess  the  impact  of  the 
ionization  rate  relaxation  time.  (The  attachment  rate  is  not  critical,  because  it  is 
not  as  sensitive  to  E-field.)  The  conclusion  reached  in  Chapter  II  is  that 
relaxation  times  below  0.1  ns  were  too  fast  for  non-equilibrium  effects  to  be 
important.  The  BVB  results  can  be  used  as  a  baseline  for  comparison.  We  need 
an  initial  electron  number  density,  Ne,  for  the  kinetic  breakdown  model. 
Following  the  BVB  analysis  of  CO 2,  we  can  assume  the  photoemitted  electrons 
from  typical  cathode  materials  are  initially  at  a  few  eV.  That  corresponds  to 
initial  velocities  (assuming  hemispherical  emission)  in  the  regime  of  10®  m/s.  If 
air  behaves  similarly  to  CO2,  we  can  expect  photoemission  currents  of  about 
10'^  A/cm2,  giving  an  initial  electron  density  of  about  10®  /cm®. 

A  preliminary  discharge  kinetic  analysis  can  be  done  with  a  nominal  external 
circuit  model.  Assuming  a  0.1— mm  gap  for  3  cm  (length)  and  0.3— cm  depth,  the 
micro-gap  "stray"  capacitance  should  be  about  10‘i‘  F  (10  pF),  and  the  stray 
inductance  is  about  3><10’i2  H.  The  complete  system  response  is  difficult  to  assess, 
but  a  total  capacitance  of  at  least  1  nF  is  reasonable  for  typical  system 
dimensions.  An  induced  voltage  of  1  kV  across  a  0.1— mm  gap  will  produce  the 
nominal  100  kV/cm  (10  MV/m)  we  need  for  fast  breakdown  to  occur. 

Using  those  external  circuit  parameters,  we  ran  a  kinetics  simulation  to  try  to 
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understand  the  likely  breakdown  behavior.  We  used  an  initial  lO^o  e/cm^  for 
convenience.  The  trend  is  easily  extrapolated  back  to  10*  or  any  other  assumed 
initial  Ne.  The  surprising  results  are  summarized  in  Figure  20.  It  shows 

essentially  a  cold  Townsend  avalanche,  with  electrons  multiplying  at  a  nearly  fixed 
rate  until  the  drawn  current  causes  voltage  collapse  (our  definition  of 
"breakdown").  That  limits  the  maximum  electron  density  to  about  10 /cm*. 
More  interestingly,  the  discharge  is  essentially  cold,  with  only  a  few  tens  of 
degrees  heating,  and  exhibits  less  complex  chemistry  than  the  transitions  to  full 
arcs  modeled  by  Rodriguez,  et  al  [Ref.  2] .  This  insight  may  lead  to  useful 
simplifications  in  a  full  micro-crack  breakdown  model. 

However,  the  breakdown  could  be  significantly  enhanced  by  planar  ionization  waves. 
In  Phase  II,  it  is  reasonable  to  build  a  1— dimensional  model  of  the  formation 
phase.  This  ID  model  would  be  based  on  solving  the  general  current  conservation 
equation 


V-CJ+f))  =  0 


where  J  =  (tE  is  conduction  current,  and  D  =  5(£E)/5t  is  displacement  current. 
In  ID,  the  gradient  operator  V-()  -♦  5()/5x.  The  model  should  include  three 
fluids  (electrons,  positive  ions  and  negative  ions)  with  their  respective  continuity 
equations: 


di 


Ne(l'i-t'a)  ••• 


di 

dN-  _ 
di 


Nel^i  •  •  • 
Nel^a  * '  • 


The  conduction  current  will  be  dominated  by  the  more  mobile  electrons,  so 
simplifications  in  the  ion  equations  are  justified.  The  general  current  conservation 
solution  could  obviate  the  need  to  solve  one  of  the  ion  continuity  equations,  since 
the  space  charge  p  is  given  by  either  V-(£E)  or  e(N+— Ng— N.).  This  has  subtle 
numerical  advantages  over  solving  the  three  continuity  equations  first,  and  then 


45 


solving  the  El— field  from  Poisson’s  equation 


V-(£E)  =  e  (N*  -  Ne  -  N.) 

since  it  avoids  differencing  errors  when  N+  a  Ne  +  N.  and  permits  numerical 
schemes  which  directly  address  the  stiff  nonlinear  coupling  between  the  E— field  and 
the  differential  conductivity  resulting  from  the  high  sensitivity  of  the  Townsend 
ionization  rate  on  E— field.  In  fact,  we  expect  the  ID  model  to  demonstrate  the 
formation  of  either  a  plane  ionization  wave  analogous  to  a  filamentary  streamer,  or 
a  non-equilibrium  Townsend  avalanche,  depending  on  the  applied  signal. 

The  energy  conservation  equation  may  not  be  strictly  necessary  in  the  QEM 
formulation,  but  it  is  a  simple  matter  to  include  it.  Electron  multiplication  terms 
of  the  form  — Ve(t'i— i^a)  or  can  easily  be  included  in  the  continuity 

and  energy  conservation  equations,  respectively.  Convective  terms  can  be 
considered  part  of  the  total  time  differential,  e.g.: 

^  =  If  +  =  (gains)  -  (losses) 

The  electron  pressure  terra  V-(NekTe)/Ne  can  be  critical  if  a  negative  glow  forms, 
where  it  dominates  over  the  E— field  force. 

Boundary  conditions  and  initial  conditions  are  straightforward.  The  anode  can  be 
modeled  as  a  perfect  electron  absorber.  The  cathode  is  dominated  by  photo- 
emission  and  ion-impact  emission  of  electrons.  The  emitted  electron  temperature 
is  not  expected  to  be  critical,  since  the  EEDF  is  quickly  dominated  by  collisions 
(in  the  first  few  /ma,  if  the  BVB  results  are  any  guide).  The  initial  conditions 
can  be  set  to  a  background  ionization  level,  and  we  let  the  voltage  rise  from  zero 
using  a  realistic  pulse  shape. 

The  envisioned  ID  model  should  yield  sufficient  insight  into  the  pre— breakdown 
mechanisms  to  explain  the  main  sources  of  variability  in  the  observed  breakdown 
behavior.  It  should  also  have  quantitative  predictive  powers  to  reproduce  the 
observable  behavior  of  the  controlled  Phase  II  experiments.  Specifically,  it  will 
model  the  cumulative  effects  between  driving-signal  half-cycles,  such  as  cumulative 
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heating  of  the  high-energy  tail  of  the  EEDF  and  the  effect  of  that  on  electron 
multiplication.  By  including  a  simplified  version  of  the  air  breakdown  kinetics,  the 
model  could  also  follow  the  cumulative  heating  leading  to  thermal  ionization. 
Thus,  a  single  code  would  model  the  formation,  heating  and  transition  to  arc 
phases,  tracking  cumulative  effects  through  many  cycles  of  the  driving  signal. 
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in  Dry  Air 


Figure  20.  Populations  in  Aperture  Breakdown  Simulation 
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VII.  PHASE  II  PROGRAM  PLAN. 


A.  PHASE  II  TECHNICAL  OBJECTIVES. 

The  overall  goal  of  the  Phase  II  Program  will  be  to  develop  a  deliverable  high- 
frequency  (HF)  air  breakdown  module  appbcable  to  micro— gaps  in  structures 
illuminated  by  electromagnetic  pulses  (EMP).  That  goal  is  supported  by  three 
objectives,  and  six  technical  tasks,  each  tied  to  a  specific  objective  as  follows: 

Objective:  Model  Development 

Task  1.  Formation  Phase  (1— D) 

Task  2.  Heating  Phase 

Objective:  Model  Validation 

Task  3.  Marginal  Effects  Study 

Task  4.  Experimental  Verification 

Objective:  Code  Integration 

Task  5.  Module  Documentation  and  Packaging 
Task  6.  Integrate  into  an  EMP  Coupling  Code 

Model  development  is  the  core  of  the  program.  The  basis  will  be  a  1— D  model 
of  a  cathode  fall  developing  into  a  Townsend  avalanche  or  ionization  wave.  The 
predominantly  local  (almost  0— D)  kinetics  describing  the  heating  phase  which  leads 
to  voltage  collapse  (i.e.  "breakdown")  will  be  built  upon  the  ID  model. 

Model  validation  assures  the  correctness  of  the  complete  model.  The  final 
criterion  is  whether  the  model  correctly  predicts  the  behavior  of  empirically 
observable  parameters,  such  as  I/V  traces  or  radiation  output.  One  approach  to 
validation  is  to  examine  neglected  effects  analytically,  which  is  the  purpose  of 
Task  3.  The  other  approach  is  to  compare  code  predictions  with  experimental 
measurements,  which  is  the  purpose  of  Task  4. 

Code  integration  assures  the  final  product  is  usable  in  a  real  application.  Task  5 
embodies  all  actions  required  to  transfer  the  module  to  a  "user".  In  Task  6,  we 
will  demonstrate  the  use  of  the  module  in  an  actual  EMP  coupling  code.  We 
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have  chosen  the  POLYANA  code  developed  by  BDM.  Thus,  this  task  requires 
participation  by  BDM  as  a  subcontractor  to  integrate  the  module  into  the 
POLYANA  code  and  test  its  performance.  BDM  will  retain  the  final  module  in 
POLYANA  for  their  use,  but  the  module  will  remain  a  separate,  marketable 
Phase  III  product. 

B.  PHASE  II  WORK  PLAN. 

Task  1.  Formation  Phase  ID  Model. 

The  technical  issues  include  the  modeling  of  cathode— fall— like  boundary  conditions, 
formation  of  a  Townsend  avalanche  and  for  a  propagating  ionization  wave,  AC 
effects,  the  effect  of  the  nonequilibrium  current  integral,  and  numerical  stiffness  of 
the  resulting  differential  equation  system.  Our  approach  is  to  extend  Tetra’s 
ELFID  code.  ELFiD  is  a  finite— volumes  electrostatic  code  which  solves  the 
generalized  current  continuity  equation: 

7-(J+D)=0 

where  J  is  conduction  current  ((tE  in  a  conductivity  model)  and  D  is  displacement 
vector  (cE).  The  finite— volumes  method  has  proven  to  be  very  forgiving  of  the 
shock— like  front  which  evolves  in  any  ionization  wave  [Ref.  21].  The  method 
captures  the  virtual  discontinuity  in  E— field,  with  no  errors  propagating  to  other 
grid  points.  ELFID  is  set  up  to  handle  planar,  cylindrical  and  spherical 

geometries,  and  could  be  modified  to  model  a  filamentary  streamer. 

The  model  should  include  three  charged  fluids  (electrons,  positive  ions  and  negative 
ions)  with  their  respective  continuity  equations.  The  conduction  current  will  be 
dominated  by  the  more  mobile  electrons;  therefore  simplifications  in  the  ion 
equations  are  justified.  The  general  current  conservation  solution  could  obviate  the 
need  to  solve  one  of  the  ion  continuity  equations,  since  the  space  charge  p  is 
given  by  either  v-(£E)  or  e(N*— Ng— N.).  This  has  subtle  numerical  advantages 
over  solving  the  three  continuity  equations  first,  and  then  solving  the  E— field  from 
Poisson’s  equation,  since  it  avoids  differencing  errors  when  N*  «  Ng+N.  and 
permits  numerical  schemes  which  directly  address  the  stiff  nonlinear  coupling 
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between  the  E— field  and  the  differential  conductivity  resulting  from  the  high 
sensitivity  of  the  Townsend  ionization  rate  on  E— field.  We  expect  the  ID  model 
to  demonstrate  the  formation  of  either  a  plane  ionization  wave  analogous  to  a 
filamentary  streamer,  or  a  nonequilibrium  Townsend  avalanche,  depending  on  the 
applied  signal. 

The  energy  conservation  equation  may  not  be  strictly  necessary  in  the  QEM 

formulation,  but  it  is  a  simple  matter  to  include  it.  Electron  multiplication  terms 
of  the  form  — Ve(i^i— i^a)  or  —<  can  easily  be  included  in  the  continuity 

and  energy  conservation  equations,  respectively.  Convective  terms  can  be 
considered  part  of  the  total  time  differential,  e.g.: 

g-j-  =  v*NV  =  (gains)  —  (losses) 

The  electron  pressure  term  7'(NeTe)/Ne  can  be  critical  if  a  negative  glow  forms, 

where  it  dominates  over  the  El— field  force. 

Boundary  conditions  and  initial  conditions  are  straightforward.  The  anode  can  be 
modeled  as  a  perfect  electron  absorber.  The  cathode  is  dominated  by 
photoemission  and  ion  impact  emission  of  electrons.  The  emitted  electron 
temperature  is  not  expected  to  be  critical,  since  the  electron  energy  distribution 

function  (EEDF)  is  quickly  dominated  by  collisions,  probably  in  the  first  few  fan. 
The  initial  conditions  can  be  set  to  a  background  ionization  level,  and  we  can  let 
the  voltage  rise  from  zero  using  a  realistic  pulse  shape. 

Task  2.  Heating  Phase  Model. 

The  Phase  II  ID  model  should  yield  sufficient  insight  into  the  pre— breakdown 

mechanisms  to  explain  the  main  sources  of  variability  in  the  observed  breakdown 

behavior.  It  should  also  have  quantitative  predictive  powers  to  reproduce  the 

observable  behavior  of  the  controlled  Phase  II  experiments.  Specifically,  it  will 
model  the  cumulative  effects  between  driving— signal  half— cycles,  such  as  cumulative 
heating  of  the  EEDF  and  the  effect  of  that  on  electron  multiplication.  By 

including  a  simplified  version  of  the  air  breakdown  kinetics,  the  model  could  also 
follow  the  cumulative  heating  leading  to  thermal  ionization.  Thus,  a  single  code 
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would  model  the  formation,  heating  and  transition  to  arc  phases,  tracking 

cumulative  effects  through  many  cycles  of  the  driving  signal. 

This  will  be  accomplished  by  adding  a  few  neutrals  to  the  kinetic  set.  However, 
since  the  neutrals  are  not  moved  by  the  E— field,  this  kinetic  system  is  local 

(0— dimensional).  The  only  spatial  influence  will  be  indirectly  through  the  electron 
fluid.  Thus,  we  will  add  a  local  model  for  few  neutrals  to  the  3— fluid  charged 
particles  model  developed  in  Task  1. 

One  issue  to  be  explored  is  methods  to  accelerate  integration  of  the  stiff  ordinary 
differential  equation  (ODE)  set.  One  alternative  is  to  change  the  independent 

variables  to  log  of  number  densities.  This  should  work  well  for  integrating 
electrons  when  they  avalanche,  but  its  effect  on  the  total  ODE  system  needs  to 

be  explored. 

Another  issue  to  be  assessed  is  the  effect  of  N2(A)  population  on  the  Boltzmann 
analysis.  As  the  NjCA)  population  increases,  superelastic  collisions  heat  the  EEDF 
in  an  energy-selective  manner.  That  increases  the  excitation  and  ionization  rates. 
This  subtle  auto-feedback  may  change  the  heating  phase  behavior  in  unexpected 
ways. 

As  these  technical  issues  are  resolved,  the  resulting  local  (0— D)  kinetics  model  will 
be  integrated  into  the  1— D  code.  Thus,  a  single  code  will  model  both  the 
formation  and  heating  phases  leading  to  breakdown. 

Task  3.  Marginal  Effects  Study. 

By  "marginal",  we  mean  effects  we  expect  to  be  negligible  to  first  order.  In  this 
task,  we  will  examine  at  least  two  such  issues:  arc  formation  and  edge  effects. 
After  breakdown,  a  high-pressure  discharge  usually  develops  into  an  arc.  On  the 
slower  time  scale  of  arc  formation  {>lfjs),  several  phenomena  ignored  in  our 
breakdown  model  become  important.  Two  phenomena  which  may  also  affect  the 
breakdown  are  pressure  buildup  and  radiation. 

In  a  constricted  (filamentary)  arc,  the  high  temperatures  cause  thermal  expansion 
which  reduces  the  neutral  number  density  N,  such  that  E/N  increases  for  a 
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constant  voltage.  In  effect,  the  critical  voltage  for  electron  avalanching  decreases. 
The  ohmic  heating  is  balanced  by  adiabatic  expansion  and  radiative  cooling  at 
equilibrium.  For  a  micro— gap,  we  don’t  expect  filamentary  streamers  to  form, 
because  the  gap  size  is  comparable  to,  or  smaller  than  the  expected  streamer 
radius.  Therefore,  we  expect  expansion  and  radiation  to  be  negligible  to  first 
order,  but  this  task  will  confirm  the  adequacy  of  that  assumption. 

Even  with  a  planar  discharge,  edge  effects  may  still  play  a  role.  The  field 
enhancements  at  the  rough  edges  of  a  micro— gap  should  cause  breakdown  to  occur 
at  the  edges  before  it  occurs  in  the  bulk  interior.  We  will  use  simple  models  to 
determine  whether  such  edge  effects  will  alter  the  overall  bulk  breakdown  times. 

These  studies  will  determine  whether  the  simple  1— D  model  is  adequate,  or  if 
thermal  expansion  or  edge  effects  must  be  included.  If  necessary,  we  will  include 
first-order  corrections  for  those  "marginal"  effects  judged  to  be  important  after  all. 

Task  4.  Experimental  Verification. 

The  ultimate  verification  of  any  computational  model  is  how  well  its  predictions 
match  experimental  observations.  To  this  end,  we  will  perform  experiments  with 
simplified  mock-ups  of  micro-gaps  in  order  to  control  the  experimental  conditions. 
The  hardest  parameter  to  control  is  the  gap  spacing,  since  in  practice  these  are 

unintentional  gaps.  Our  experimental  approach  is  predicated  on  pressure  scaling. 
By  conducting  the  experiments  at  low  pressures,  we  can  scale  all  dimensions  larger 
in  inverse  proportion  to  the  pressure.  For  example,  by  conducting  the  experiments 
at  10 '2  atm  (7.6  Torr),  we  can  build  a  gap  mock-up  100  times  larger  than  what 
is  being  modeled.  As  a  side  benefit,  all  time  scales  are  slowed  down  by  the  same 
factor,  relieving  the  data  acquisition  requirements. 

Two  laboratory  experiments  will  be  designed  and  fabricated  to  collect  the  data. 

One  experiment  utilizes  a  transmission  line  approach  and  the  other  a  lumped 

circuit  approach.  The  difference  is  that  the  transmission  line  will  propagate  a 
double— exponential  electromagnetic  field  across  the  slot  that  will  induce  electric  and 
magnetic  dipoles  at  the  slot.  In  the  lumped  circuit,  the  double-exponential 
electric  field  will  be  applied  directly  across  the  slot. 
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The  transmission  line  approach  is  shown  in  Figure  21.  A  double-exponential 
voltage  is  produced  and  propagated  down  a  transmission  line  to  a  matched  load. 
A  mismatched  system  can  be  employed  to  determine  how  oscillating  currents  affect 
the  discharge  characteristics  of  the  slot.  A  slot  is  in  the  ground  plane  portion  of 
the  transmission  line.  The  slot  is  thin  and  nearly  spans  the  entire  width  of  the 
transmission  line.  Data  of  interest  include  voltage  across  the  slot,  current  through 
the  slot,  open— shutter  photographs  of  the  discharge,  and  time  resolved  UV 
measurements.  The  voltage  and  current  data  will  provide  the  characteristics  of 
voltage  collapse  and  current  growth,  respectively.  The  closure  time  (time  for  the 
air  to  become  electrically  conductive)  for  the  slot  will  be  determined  from  the 
voltage  and  current  waveform.  Open-Shutter  photographs  may  provide  information 
on  whether  the  discharge  is  filamentary  or  diffuse,  which  affects  the  breakdovm 
mechanisms  involved.  A  temporal  profile  of  the  UV  radiation  emitted  by  the 
discharge  will  be  measured.  From  this  measurement,  the  electron  number  density 
can  be  determined.  (The  electron  number  density  is  an  important  quantity  in  the 
theoretical  analysis.)  Required  equipment  include  high  frequency  current  and 
voltage  probes,  a  500  MS/s  digital  storage  oscilloscope  (three  channels  preferred), 
open-shutter  camera,  UV  detector  circuit,  power  supply,  circuit  components 
(capacitors,  inductors  and  resistors),  high  vacuum  contact  switch,  copper  sheet,  and 
other  miscellaneous  materials. 

The  lumped— circuit  approach  is  shown  in  Figure  22.  A  double— exponential  voltage 
is  simply  produced  across  the  slot.  The  data  of  interest  and  required  equipment 
are  the  same  as  for  the  transmission-line  approach  except  that  the  circuit 
components  are  different. 

Experiments  will  be  conducted  down  to  about  10  Torr.  Experimental  variables 
and  parameters  include  relative  humidity,  type  of  air  (natural  or  artificial),  and 
pressure.  All  materials  in  the  vessel  will  be  waterproof  and  rust— resistant. 
Required  equipment  includes  an  airtight  vessel,  roughing  pump,  humidity  meter, 
thermometer,  and  other  miscellaneous  materials. 

In  designing  the  experiments,  we  will  use  computer  simulations  of  the  experimental 
circuits  with  an  open  circuit  gap  to  determine  the  applied  signals.  We  will 
produce  mechanical  drawings  of  the  transmission  line  and  the  airtight  vessel.  The 
type  and  requirements  for  the  sensors  will  be  determined  from  results  of  the 
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drcmt  simulations.  A  test  plan  will  be  developed  that  details  the  type  and 
volume  of  data  needed  to  support  the  computational  analysis. 

After  ordering  the  parts  and  materials,  we  will  fabricate  and  check  out  the 
experiment.  We  will  build  the  experiment  and  verify  operation  of  the  sensors  and 
recording  instrumentation.  The  collected  experimental  data  will  include  the 
experimental  pairameters  of  pressure,  relative  humidity,  type  of  air,  and  electric 
field  across  the  slot. 

We  will  analyze  the  experimental  data  and  compare  it  with  computational  results. 
Here  our  analysis  will  include  correlation  of  discharge  conduction  time  with  spatial 
and  temporal  profiles  of  the  discharge. 

Numerical  simulations  of  the  experimental  conditions  will  be  performed  for  direct 
comparison  with  experimental  measurements.  The  simulations  will  include  external 
circuit  models  matching  the  experimental  driving  circuits.  Output  parameters  of 
interest  include  the  current  and  voltage  traces,  closing  times,  and  photodetector 
signal  compared  to  computed  excitation  rates.  Favorable  comparisons  will  then 
establish  confidence  in  the  computational  model. 

Task  5.  Module  Documentation  and  Packaging. 

In  order  to  permit  a  third  party  user  to  use  the  resulting  module,  it  is  imperative 
that  it  be  packaged  in  a  user-friendly  form.  This  includes  documentation  of  the 
basic  module  itself  and  auxiliary  programs. 

During  the  development  of  the  code,  we  will  develop  drivers  for  the  module  that 
exercises  it  (and  perhaps  sub-modules).  In  addition,  we  will  also  develop  post¬ 
processors  that  read,  analyze  and  display  diagnostic  data  pertaining  to  the  status 
of  the  module’s  internal  data  structure.  Some  of  these  drivers  and  post— processors 
will  be  useful  to  a  user  of  the  module,  and  need  to  be  packaged  accordingly. 

Most  importantly,  a  User  Manual  will  be  written  for  the  basic  module  and  for  the 
drivers  and  post— processors.  It  will  contain  detailed  instructions  on  how  to  use 
them,  down  to  the  level  of  defining  each  argument  in  the  call  sequence  in  the 
basic  module.  This  is  indispensable  for  a  programmer  to  integrate  the  module 
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into  a  larger  EM  coupling  code. 


Task  6.  Integrate  Into  An  EMP  Coupling  Code. 

The  ultimate  test  of  whether  the  Phase  II  module  is  usable  in  an  EMP  coupling 
code  is  to  integrate  it  into  one.  Otherwise,  there  may  be  code  coupling  problems 
not  obvious  a  priori  which  would  surface  only  when  the  first  real  user  tries  to 
use  the  module.  For  this  purpose,  we  have  chosen  the  POLYANA  code  developed 
by  BDM  International,  Inc. 

BDM’s  POLYANA  code  offers  a  framework  which  is  flexible,  and  works  in  the  time 
domain.  First,  we  will  give  a  brief  description  of  POLYANA  then  raise  some  issues 
to  be  addressed.  Finally,  we  present  an  example  of  how  we  might  address 
integration  of  the  QEM  module  and  data  (I/O)  handoffs.  POLYANA  is  a  time- 
dependent  (explicit  time— domain)  finite  difference  code  for  solving  Maxwell’s  curl 
equations  in  a  three  dimensional  geometry.  The  code  is  capable  of  modeling 
complex  objects  by  specifying  the  scalar  conductivity  and  dielectric  constant  at 
each  point  in  the  problem  space.  Small  sub— grid  features  such  as  thin  wires  and 
struts  can,  and  have  been,  nodeled  using  a  built-in  line  integral  formalism. 

POLYANA  also  has  an  optional  air  chemistry  model  for  low  altitude  source  region 
coupling.  Scattered  radiation  is  terminated  at  the  edge  of  the  problem  space  by  a 
radiation  absorbing  boundary  condition.  The  user  specifies  the  finite  difference 
grid,  material  properties,  surfaces  for  the  Huygen’s  sources,  and  sub— grid  structure 
locations.  It  is  written  in  highly  transportable  FORTRAN  77  which  can  run  on  a 
CRAY,  VAX,  Sun  4/110,  or  386/486  class  PC  with  an  appropriate  amount  of 
RAM.  POLYANA  is  a  BDM  constructed  intellectual  descendant  of  the  old  THRED 
series  codes  and  has  itself  been  the  subject  of  several  upgrades  and/or  special 

purpose  modifications.  It  has  been  used  in  the  past  primarily  for  calculating  the 
EMP  coupling  to  a  wide  variety  of  structures  including  the  Peacekeepe’  missile 
and  the  Rail  Garrison  system. 

Initially,  we  propose  to  use  a  cartesian  coordinate  version  of  POLYANA  rather  than 

a  conformal  one  for  our  initial  effort.  We  can  change  to  numerically  generated 

coordinates  later.  The  initial  concern  is  timing,  i.e.  how  fast  this  modified  code 
will  run  a  given  problem  with  the  QEM  module  added.  Several  issues  need  to  be 
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resolved  in  order  to  determine  how  to  proceed  in  detail,  e.g.: 


Do  we  break  the  problem  into  three  pieces  and  use  prescribed  sources  to 
connect  the  pieces? 

If  so,  should  the  prescribed  sources  be  consistent? 

If  we  don’t  break  the  problem  into  pieces  how  do  we  want  to  handle  the 
sub— grid  time  resolution  problem  in  the  POLYANA  interface? 

Specifically,  how  do  we  address  the  stochastic  part  of  the  problem  relative 
to  the  deterministic  field  calculations? 

Our  current  thinking  on  these  issues  is  as  follows:  Self-consistency  is  highly 
desirable,  and  should  be  maintained,  unless  it  makes  the  coupling  burdensome.  In 
a  stand-alone  mode,  the  QEM  module  would  be  coupled  to  a  simple  external 
circuit  model,  so  currents  and  voltages  can  be  solved  self— consistently.  Ideally,  the 
coupled  module  would  be  structured  so  that  the  external  code  can  model  a 
complete  aperture  as  a  coupled  network  of  micro— gaps  or  cells,  each  treated  as  a 
1— D  discharge  independent  of  all  the  others  except  through  the  external  circuit’s 
coupling. 

One  simple  linearization  scheme  is  to  let  the  external  code  prescribe  voltage  during 
a  time  step  of  its  choosing.  The  QEM  module  would  then  compute  current  (at 
the  end  of  the  time  step  and  integrated  over  the  time  step),  using  internal  time 
sub-steps  if  needed.  That  way  the  external  code  time-step  management  is  totally 
independent  of  the  QEM  internal  time  management.  In  practice,  the  need  for 
CPU  expenditure  control  will  dictate  some  compromise  or  optional  choices,  perhaps 
patterned  after  the  Livermore  Solver  of  Ordinary  Differential  Equations  (LSODE) 
time  controls  [Ref.  22]. 

Higher  order  linearization  schemes  are  possible,  based  on  a  simple  external  circuit 
model  "local"  to  each  micro-gap  cell,  tied  together  by  the  external  code.  The 
nonlinear  nature  of  the  problem  dictates  that  self-consistency  be  handled  in  an 
outer  loop  under  the  external  code’s  control.  A  robust  linearization  scheme  will 
assure  fast  convergence  in  the  regime  where  the  micro— gap  is  nearly  linear 
(impedance  not  changing  rapidly),  and  yet  handle  the  fast  voltage  collapse  and 
current  rise  which  characterizes  breakdown. 
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The  stochastic  part  of  the  problem  can  be  handled  by  random  or  pseudo— random 
assignments  of  effective  gap  length  to  the  various  QEM  cells.  Each  cell  will 
represent  a  small  area  of  the  aperture  perimeter,  over  which  the  gap  can  be 
treated  as  constant.  Thus,  the  total  model  will  be  quasi— 2D  or  quasi— 3D.  The 
other  stochastic  part  is  the  initial  conditions  in  each  gap— cell.  Methods  for 
randomizing  those  will  be  evaluated  with  respect  to  physical  correctness,  flexibility, 
ease  of  use,  etc. 

The  purpose  of  this  proposed  effort  is  to  understand  how  electromagnetic  fields  go 
from  the  outside  of  a  structure  to  the  inside  through  micro— gaps  which  break 
down.  The  incident  EMP  field  is  generated  in  the  problem  space  outside  the 
enclosure  by  a  user— specified  Huygen’s  source;  then  the  self-consistent  electro¬ 

magnetic  object  response  is  computed. 

Given  that  we  have  self-consistent  fields  at  each  grid  cell  in  the  problem  space, 
let’s  then  examine  how  a  user  would  view  the  micro— gap.  A  basic  assumption  is 
that  the  gaps  are  electrically  short  since  we  are  working  a  variant  of  the  "classic" 
EMP  problem  (even  if  we  use  2169  waveforms).  After  a  random  number  draw  is 
used  to  select  which  material  grid  cell,  or  cells,  is  to  have  a  micro-gap,  the  local 
self-consistent  electric  and  magnetic  fields  (or  surface  current  and  charge  density) 
will  be  passed  to  the  QEM  subroutine  to  do  the  gap  discharge  calculations. 

POLYANA  could  pass  those  quantities  to  the  QEM  subroutine  at  each  time— step. 
Then  those  fields  could  be  modified  to  reflect  the  gap  geometry  by  an 

approximate  calculation;  e.g.  by  using  normalized  De  Meulenare— Von  Bladel 

polarizabilities  to  calculate  dipole  moments  and  then  the  modified  field  vector 
components  (or  surface  current  density  and  charge  values)  prior  to  beginning,  or 
continuing,  the  discharge  calculation. 

Consider  the  QEM  handoff  back  to  polyana.  The  QEM  subroutine  should  pass 
back  gap  currents  and  voltages  at  each  POLYANA  time-step.  The  sub-grid  line 
integral  formalism  will  reconnect  these  quantities  to  the  field  nodes  on  the  finite 
difference  grid  inside  the  macroscopic  structure.  We  then  just  need  to  tell 
POLYANA  where  to  locate  our  test  points  to  obtain  field  information  for  further 

post  processing  and  presentation. 

Clearly,  there  are  many  issues  to  be  resolved  before  the  module  interface  can  be 
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designed  in  detail.  Actually  integrating  th  QEM  module  into  a  real  EMP 
coupling  code  like  POLYANA  insures  that  the  interface  issues  are  examined  from 
both  sides. 

Finally,  several  test  cases  will  be  run  to  demonstrate  the  usefulness  of  a  nonlinear 
micro— gap  breakdown  option  in  a  full  EM  coupling  code.  The  test  cases  will  be 
defined  jointly  by  Tetra  and  BDM  in  consultation  with  the  Contract  Technical 
Monitor,  in  order  to  insure  that  they  are  feasible,  meaningful  and  relevant  to  the 
mission  of  Harry  Diamond  Labs. 

Task  7.  Management  and  Reporting. 

Progress  Reports  will  be  submitted  quarterly.  A  written  Final  Technical  Report 
will  document  model  derivations,  results  of  test  cases,  and  comparisons  of 
experimental  results  with  simulations.  The  final  QEM  module,  drivers  and  post¬ 
processors,  along  with  the  Users  Manual  will  be  considered  contract  deliverables. 
One  or  two  program  reviews  are  recommended.  A  mid-term  and  a  final 

presentation  at  Harry  Diamond  Labs  is  recommended. 

Schedule  (months) 


1.  Formation  Phase 

2.  Heating  Phase 

3.  Marginal  Effects 

4.  Experimental  Verification 

5.  Documentation  k  Packaging 

6.  Integrate  in  POLYANA 

7.  Management  k  Reporting 


0  3  6  9  12  15  18 
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Optional  tasks  could  be  defined  and  included  in  the  Phase  II  project  as  pre-costed 
but  unfunded  options  to  the  basic  contract.  One  possibility  is  to  extend  Task  3 
into  development  of  a  full  2— D  breakdown  model  for  $50K  to  $70K.  Another 
logical  extension  of  Task  3  would  be  to  develop  an  arc  phase  model  capable  of 
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predicting  recovery  times  (which  relate  to  repetition  rates).  That  may  cost  $50  to 
$60K.  A  logical  extension  of  Task  5  may  be  to  develop  user-4nendly  driving 
programs  using  Graphic  Input  (GIN)  and  other  advanced  techniques  as  user 
training  aids.  That  would  cost  between  S25K  and  $35K.  A  logical  extension  of 
Task  6  is  to  apply  the  integrated  POLYANA/QEM  code  to  analyze  a  specific 
problem  of  vital  interest  to  HDL.  Depending  on  the  difficulty  of  the  specified 
problem,  that  may  cost  between  $25K  and  $50K. 
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Figtire  21.  Transmission  line  experimental  setup. 
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Figure  22.  Lumped  circuit  experimental  setup. 
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VIII. 


SUMMARY  AND  CONCLUSIONS. 


In  Phase  I  of  the  program,  we  developed  the  basic  tools  needed  to  build  a  high- 
frequency  air  breakdown  model.  We  developed  a  Quasi— Equilibrium  Methodology 
(QEM)  for  extending  fluid  electron  transport  models  to  the  nonequilibrium  regime. 
Our  analysis  of  nonequilibrium  electron  transport  discovered  and  explained  a 
velocity  overshoot  effect  that  is  not  generally  known.  That  tendency  of  the 
electrons  to  overshoot  tne  equilibrium  drift  velocity  was  shown  to  have  a  profound 
effect  on  streamer  behavior.  We  examined  the  extent  of  possible  simplifications  to 
the  basic  air  breakdown  model,  concluding  that  factors  of  2  to  5  speed-up  are 

possible.  We  have  identified  the  major  phenomenological  issues  in  a  modeling 
strategy,  and  examined  nonequilibrium  effects  on  ionization  waves.  Finally,  we 
presented  a  plan  to  complete  development  of  a  computational  model  of  High 
Frequency  air  breakdown  suitable  for  incorporation  into  EMP  coupling  codes. 

Although  the  main  application  is  for  including  nonlinear  breakdown  phenomena  in 
EM  coupling  problems,  there  are  other  applications  for  a  High  Frequency 

breakdown  model.  The  nonequilibrium  effects  on  ionization  waves  can  be  very 
important.  An  accurate  ionization  wave  model  will  find  applications  in  arc  lamp 
and  spark  gap  design  [Ref.  16],  the  study  of  HV  vacuum  insulator  surface 

flashover  [Ref.  23],  and  perhaps  in  advanced  concepts  using  controlled  ionization 
waves  to  accelerate  non— relativistic  ion  beams  [Ref.  24].  It  also  seems  to  have 
applications  in  modeling  of  air  avalanche  switches  for  EMP  generation  [Ref.  25], 

and  the  velocity  overshoot  described  in  Chapter  II  has  been  observed  in  solid  state 
devices  [Ref.  26].  The  QEM  applied  to  other  gases  may  also  find  application  in 

low  pressure  glow  discharge  switches,  such  as  the  Cs-vapor  tadtron  [Ref.  27]. 

However,  in  order  to  limit  the  scope,  the  focus  of  the  Phase  II  effort  will  be  to 
model  micro— gap  breakdown  in  apertures  of  shielding  structures  illuminated  by 
EMP. 

The  final  product  will  be  a  computational  module  incorporating  the  air  breakdown 
model,  which  can  be  either  run  in  a  stand-alone  configuration,  or  incorporated 

into  a  larger  EMP  coupling  code.  To  demonstrate  the  usefulness  of  the  module, 

we  will  indeed  integrate  it  into  the  POLYANA  EMP  coupling  code  developed  by 
BDM.  Tetra  plans  to  market  the  module  separately  in  Phase  III,  not  only  to  the 
EMP  community,  but  also  to  the  community  designing  and  developing  glow 
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discharge  switches,  lasers  and  arc  lamps,  as  we  have  successfully  done  with  our 

Electric  Field  Analysis  (tetraELF  )  codes  in  which  the  complete  range  of 
marketing  activities  were  conducted.  This  included  market  research,  market 
analysis,  development  of  descriptive  materials  and  specification  sheets,  all  followed 
by  an  international  sales  campaign. 

As  a  direct  result  of  the  Phase  I  accomplishments  of  the  Quasi— Equilibrium  Model 
of  HF  Air  Breakdown  project,  we  are  in  a  position  to  develop  in  Phase  11  a 
breakdown  model  to  predict  breakdown  in  apertures  illuminated  by  EMP.  Such  a 
computational  model,  once  fully  verified,  would  have  applications  in  arc  lamp  and 
spark  gap  design,  including  switches  for  EMP  generation,  vacuum  HV  insulators, 
solid-state  devices,  and  probably  others  we  have  not  discovered. 

For  approximately  $300K,  HDL  would  obtain  what  very  weU  could  become  the 
standard  computational  model  for  assessing  the  effect  of  air  breakdown  on  the 
response  of  electronic  systems  to  EMP  and  other  high— power  electromagnetic 
threats.  Our  plans  to  commercially  market  the  computational  module  in  Phase  III 
would  make  this  modeling  capability  accessible  to  the  EMP  community  and  those 
in  other  application  areas. 
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Appendix  A 

Relaxation  to  a  Sliding  Quasi-Equilibrium  —  An  Algorithm 

The  purpose  of  this  Appendix  is  to  document  the  derivation  of  an  algorithm  for 
integrating  the  elementary  differential  equations  on  which  the  Quasi— Equilibrium 
Methodology  (QEM)  is  based.  The  basis  of  the  QEM  model  is  that  the  collision 
frequencies  and  rates  relax  toward  values  in  equilibrium  with  the  local  electric  field 
not  instantaneously,  but  with  certain  characteristic  lag  times. 

Consider  electron  velocity  as  an  illustration.  The  momentum  conservation  equation 
without  convective  terms  reduces  to  an  electron  fluid  acceleration  equation: 

V  =  -^  E  -  V  1/ 
m 

(V  is  electron  fluid  velocity,  e/m  is  charge  to  mass  ratio,  and  E  is  the  driving 
dectric  field.)  Under  Locd  Field  Approximation  (LFA)  methods,  the  collision 
frequency  v  driving  the  friction-dike  loss  term  is  considered  a  function  of  local, 
instantaneous  E— field  (usually  pre— tabulated  from  empirical  measurements  or 
Boltzmann  solutions).  The  QEM  approach  varies  in  that  v  is  calculated  as  the 
solution  of  an  accompanying  Ordinary  Differential  Equation  (ODE): 

I  =  (i^i/q)/r 

where  the  equilibrium  solution  Uq  and  characteristic  time  r  are  considered  functions 
of  instantaneous  E. 

Electron  number  density,  with  a  source  term  S  and  attachment  loss  rate  a,  is 
governed  by  an  ODE  of  similar  form: 

Ne  =  S  —  Ne  a 

We  may  generalize  the  ODE  form  by  calling  the  dependent  variable  Y,  the 
driving  function  S,  and  the  loss/collision  rate  R,  thus: 

Y  =  S-  YR  and  fl  =  (R-Rq)/r 

where  Rq  is  the  equilibrium  value  R  tends  toward  with  a  relaxation  time  r. 

Consider  a  discrete  time  step  At  such  that  the  driving  functions  S,  Rq  and  r  are 
well-behaved  enough  to  use  linear  interpolation.  If  Rq  and  r  were  constant,  the 
second  equation  would  have  an  exact  solution 

R(t)  =  Rq  (Ro-Rq)  exp(-t/r) 

where  Rq  h  R(0)  is  the  initial  value.  This  is  the  essential  relaxation  behavior: 
the  instantaneous  difference  from  the  equilibrium  solution  decays  exponentially  in 
time. 

A  graphical  interpretation  of  the  ODE  is  that  the  slope  of  the  R  vs  t  curve  must 
point  at  each  instant  toward  the  tip  of  an  arrow  at  the  equilibrium  value  Rq 
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delayed  by  one  time  constant  r.  If  both  are  changing  in  time,  the  system  is 
"chasing”  a  moving  target.  A  linear  interpolation  of  the  "target"  represents  a 
"Sliding  Quasi— Equilibrium",  such  as  the  long— dashed  line  in  Figure  1  labelled 
"SQE".  Decay  of  the  initi^  difference  from  that  SQE  represents  a  "Relaxation  to 
a  Sliding  Quasi— Equilibrium"  (RSQE).  If  the  initial  and  final  equilibrium 
R— values  are  Ro  and  Ri,  then  the  targets  at  0  and  At  are  easily  estimated  by 
linear  int^olation.  Let’s  designate  them  Ro’  and  Ri’.  Thus,  the  RSQE 
algorithm  is: 


Rn  =  R(At)  =  R,’  +  (Ro-RoO 

Since  we  have  chosen  At  such  that  r  is  well-behaved,  the  integral  can  be 
evaluated  by  linear  averaging: 

/(1/T)(it  »  <At/T> 

where  the  angular  brackets  <  >  stand  for  the  operation  of  taking  the  linear 
average  over  At. 

In  Figure  1,  the  function  R(t)  was  calculated  at  32  sub— steps  of  At,  and  the 
more  accurate  32-step  solution  is  plotted  along  with  the  one— step  solution  curve. 
Notice  they  are  nearly  indistinguishable. 

The  same  procedure  could  then  be  applied  to  the  primary  ODE,  with  the 
initial/final  values: 


Yo  =  So/Ro  and  Yi  =  Si/Rn 
However,  the  driving  function  R(t)  is  not  linear.  The  integral 

/  R  dt  =  <Rs>  At  +  (Ro-Ro’)  r 

where  <Rs>  =  0.5(Ro’+Ri’)  is  the  linear  average  SQE  R— value,  is  accurate  to 

the  same  order  as  R(At).  Using  linear  interpolation  to  define  SQE  target  values 
Yo’  and  Yj’  as  illustrated  in  Figure  2,  we  would  calculate  the  solution  as: 

Y(At)  s  Yi’  +  (Yo-Yo’) 

Notice  the  one— step  function  (the  curve  labeled  "RSQE"  in  Figure  2)  overestimates 
Y(t)  compared  to  the  more  accurate  32-step  solution.  The  reason  is  that  the 
true  SQE  target  depends  on  an  R(t)  function  which  is  not  linear.  By  the  time 
we  get  near  At,  R  has  relaxed  to  much  nearer  the  Rj  equilibrium  value  than 

implied  by  a  linear  interpolation.  It  is  important  that  the  FISQE  solution  not  be 

less  accurate  than  an  equilibrium  solution  for  large  At  >>  1/R. 

A  next-order  correction  is  needed  to  account  for  the  nonlinear  shape  of  the 
Y-function  SQE.  For  a  short  characteristic  time  1/R  <<  At,  it  is  reasonable 
to  back  up  to  t’  =  At— 1/R,  with  the  expectation  that  the  forward  projection 
t’4-l/Rft’)  will  be  close  to  At.  and  thus  a  better  basis  for  estimating  Yi’  at  At. 
Since  tne  analytic  function  R(t)  is  known,  we  can  calculate 

Y’(t’)  =  S(t’)/R{t’) 
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where  S(t’)  is  calculated  by  linear  interpolation.  This  is  the  basis  for  a  nonlinear 

SQE  Y-Tunction  target.  However,  for  stability,  we  must  use  a  form  for  t’  which 

does  not  result  in  a  backward  projection  (t’<  0).  Defining  Rave  =  /Rdt/At,  the 
form  t’  =  At  exp^— 1//Rdt)  tends  to  0  for  1/Rave  >>  and  tends 

asymptotically  to  At— 1/Rave  for  1/Rave  <<  At.  That  is  the  form  used  in  the 
computations  displayed  in  Fi^e  3.  Notice  that  the  modified  RSQE  one-step 
solution  is  indistinguishable  nom  the  more  accurate  32-step  solution  near  the 
beginning  and  near  the  end  of  the  test  case  time  step,  but  an  error  is  evident 

in-between.  (Keep  in  mind  that  the  test  case  is  intentionally  stressful.) 

The  source  of  the  error  is  that  the  simple  linear  SQE  was  being  used  for  the 

extrapolation  to  0  of  Yo’.  For  consistency,  we  ought  to  use  the  same  Y’,Yi  data 
pair  to  re— compute  Yq’.  That  amounts  to  using  a  near— implicit  tangent  to  the 
non^near  SQE  as  the  basis  for  the  Y— relaxation  algorithm.  It  makes  the 
algorithm  more  implicit,  by  weighting  data  near  the  end  of  the  time  step  heavier 
as  the  characteristic  time  1/Rave  becomes  small  compared  to  At.  With  that 
modification,  the  "implicit''  R5QE  algorithm  gave  the  results  shown  in  Figure  4. 
Notice  the  agreement  is  excellent  at  ^  times. 


However,  that  does  not  guarantee  the  1— step  integral  of  Y  is  acceptably  accurate. 
Within  the  RSQE  formalism,  the  integral  is  calculated  as: 


/  Y  dt  «  At 


<Ys>  +  AYo 


where  <Ys>  is  the  average  of  the  intercepts  Yo’  and  Yj’,  AYo  is  the  initial 
difference  to  the  intercept  (Yq- Yo’).  The  problem  with  that  algorithm  is 

illustrated  well  by  a  test  case  designed  to  accentuate  the  tendency  of  Y(t)  to 
overshoot  its  eventual  equilibrium  value.  The  test  case  uses  initial  values  Yo  = 
0,  and  Ro  =  1,  tending  toward  a  constant  equilibrium  value  Rq  =  5  in  a 
constant  relaxation  time  r  =  10,  with  a  constant  source  term  S  =  5.  The 
solution  Y  vs  t  curve  is  shown  in  Figure  5. 


The  simple  "linear"  RSQE  algorithm  seriously  over-estimates  the  integral,  but  the 
non-linear  algorithm  using  the  "implicit"  offset  t’  near  At— 1/Rave  under-estimates 
it.  The  reason  is  that  a  tangent  to  the  nonlinear  SQE  near  the  end  of  the  time 
step  extrapolated  back  to  t  =  0  misses  the  peak.  An  alternate  algorithm  using  a 
more  "explicit"  displaced  time  t’  near  I/Rq  over-estimates  the  integral,  but  not  as 
badly  as  the  linear  algorithm.  A  new  approach  is  needed. 


We  decided  to  estimate  the  "dominant"  time  tj  at  which  we  should  expect 

Y(td)  At  «  /  Y  dt 


Substituting  the  RSQE  algorithm  expressions,  we  get 


Y,  +  AY.  »  <Y.>  +  AY. 

where  Yg  and  <Ys>  are  respectively  the  instantaneous  and  average  SQE  values. 
Assuming  Ys(td)  «  <Ys>, 
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where  we  used  the  e^lidt  time  constant  1/Ro  outside  the  bracket,  to  insure  that 
the  peak  region  is  given  adequate  weight.  As  shown  in  Figure  6,  the  resulting 
solution  for  the  test  case  is  nearly  indistinguishable  from  the  "correct"  32-step 
calculation  (dashed  line).  Notice  all  the  algorithms  are  good  at  small  times,  but 
the  rejected  algorithms  have  errors  which  grow  with  time. 

In  conclusion,  an  algorithm  has  been  developed  to  integrate  the  ODE  pair 

t  =  S  -  Y  R  and  R  =  (R-Rq)/r 

by  a  method  we  call  Relaxation  to  a  Sliding  Quasi— EquiUbrium.  It  has  the 
correct  asymptotic  behavior  for  time  steps  At  much  smaller  than  or  much  greater 
than  the  dominant  characteristic  times,  and  behaves  well  anywhere  in-between. 
The  ODE  characteristic  times  r  and  1/R  may  well  be  much  smaller  than  the 
characteristic  time  At  on  which  the  driving  functions  S,  Rq  and  r  change.  The 
RSQE  algorithm  applies  even  for  At  >>  1/R,  a  situation  which  is  too  stiff  for 
conventional  numerical  integration  methods.  (For  example,  if  R  At  «  10^,  a 
"brute  force"  algorithm  would  need  10^  substeps  for  a  1%  tolerance.) 
Furthermore,  RSQE  time  steps  can  be  chosen  by  the  more  bberal  criterion  of 
limiting  changes  in  the  slopes  of  driving  functions,  rather  than  the  more  restrictive 
criterion  of  bmiting  changes  in  the  functions  themselves. 

In  the  long  rua,  the  RSQE  algorithm  is  not  merely  more  efficient;  it  may  be  the 
only  feasible  method  to  handle  the  strong  stiffness  introduced  by  the 
Quasi-Equilibr  im  Methodology. 
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Industrial  environments,  at  least  outside  of  clean  rooms,  are  contaminated  by 
airborne  particulates  and  particulates  on  most  surfaces.  Airborne  particulates, 
whether  they  are  dielectrics  or  conductors,  are  capable  of  acquiring  and  sustaining 
an  electrical  charge.  In  this  respect,  particulates  can  often  be  active  components 
of  plasmas  sustained  in  ambient  air. 

Large  ungrounded  particles  in  plasmas  (sizes  commensurate  to  the  Debye  length) 
will  most  often  negatively  charge  in  the  same  fashion  that  a  floating  surface  in 
contact  with  a  plasma  will  negatively  charge  to  balance  the  flux  of  electrons  and 
ions  to  its  surface.  In  this  respect,  particulates  resemble  massively  large,  multiply 
charge  negative  ions.  Particles  having  sizes  of  a  few  microns  will  have  100s  to 
1000s  of  elementary  charges  on  their  surfaces. 

The  fact  that  particles  are  suspended  in  the  plasma  (either  electrostatically  or 
buoyantly)  and  are  negatively  charged  detrimentally  impacts  the  plasma.  The 
consequences  of  particulate  contamination  of  plasmas  are  at  least  four-fold.  First, 
the  electron  energy  distribution  (EED)  of  a  contaminated  plasma  is  shifted  to 
lower  energies  compared  to  the  EED  in  a  pristine  plasma  at  the  same  E/N 
(electric  field/number  density)  (Ref.  1).  Second,  due  to  their  effective  large 

momentum  transfer  cross  sections,  particulates  will  channel  current  in  the  plasma 
out  of  regions  of  high  contamination  into  regions  of  low  contamination  (Ref.  2). 
Third,  particles  are  large  recombination  centers  for  electrons  and  ions  which  can 
significantly  affect  the  balance  between  electron  production  and  electron  loss. 
Fourth,  during  transient  operation,  the  number  of  electrons  required  to  charge  a 
particle  may  be  a  non— neglible  fraction  of  the  total  charge  available. 

To  determine  whether  particles  must  be  considered  during  the  development  and 
propagation  of  streamers,  the  following  scaling  laws  should  be  considered. 

1.  Assume  that  the  streamer  is  propagating  into  an  otherwise  nonionized  gas 
which  is  contaminated  by  particles  having  density  N  and  radius  R.  Further 
assume  that  the  path  length  of  propagation  is  1.  A  streamer  will  certainly 
encounter  these  uncharged  particles  if  the  fractional  area  density  exceeds  unity. 
That  is,  particulate  contamination  is  important  for  these  conditions  if 

47rR2Nf  >  1. 

2.  If  the  streamer  is  propagating  into  an  ionized  plasma,  the  critical  radius 
increases  to  approximately  R+A,  where  A  is  the  Debye  length.  [A  =  750  * 

(Tg/ne)^  cm,  Tg  =  electron  temperature  (eV),  Ug  =  electron  density  (cm'3)]. 
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3.  The  charge  required  to  raise  the  potential  of  a  spherical  dielectric  having 
radius  R  to  V  volts  is  V*4Teo*R.-  The  rate  of  charging  is  approximately  j-A, 
where  j  is  the  current  density  and  A  is  the  cross  sectional  area  of  the  particle. 
Assuming  the  particle  will  charge  to  a  few  times  the  electron  temperature,  the 
charging  time  is 


At  =  2.5*  eo’Te/(j*R) 

4.  The  effective  rate  of  recombination  of  electrons  and  ions  on  the  surface  of  a 
particle  is  given  by  the  flux  of  ions  to  its  surface,  since  all  ions  will  neutralize 
once  they  strike  the  negative  charge.  The  effective  rate  coefficient  for 
recombination  is  therefore  given  by  the  cross  sectional  area  of  the  particle  and  the 
ion  thermal  velocity. 


kr  =  Vion(^R^)  cm^s'i 

The  total  rate  of  recombination  is  Vr  =  kr-N,  and  can  exceed  lO^  s'l  in 
moderately  hot  and  contaminated  plasmas. 

(NOTE:  These  expressions  must  be  evaluated  for  the  conditions  of  a  particular 

problem.) 

In  previous  studies  of  the  effects  of  particles  in  low  pressure  glow  discharges,  it 
was  determined  that  the  operating  E/N  of  the  discharge  increased  with  increasing 
particle  contamination  because  the  rate  coefficients  for  ionization  decreased  (Ref.  1 
and  2).  In  the  0.1  -  10  Torr  range,  significant  effects  were  found  for  particulate 
densities  (radius  1  /xra)  exceeding  10^  cm*^.  In  general,  the  threshold  for  this 
effect  increased  with  increasing  pressure  and  increasing  E/N.  As  a  practical 
consideration,  in  inhomogeneously  contaminated  plasmas  this  may  not  be  a 
particularly  important  effect.  This  apparent  contradiction  results  from  the  fact 
that  once  the  particles  charge,  the  current  will  be  channeled  away  from  heavily 
contaminated  regions  into  less  contaminated  regions.  Hence,  the  local  lowering  of 
the  ionization  rate  coefficients  may  not  be  particularly  important.  In  high 
pressure  plasmas  where  the  local  field  approximation  is  v^d  (p  >  10s  — 
100  Torr)  and  the  Debye  length  is  short  (/zm’s)  the  isolation  of  the  particles  by 
these  mechanisms  is  fairly  efficient. 

The  fact  that  current  is  constricted  into  less  contaminated  regions  will  have 
secondaiy  effects.  For  example,  in  electronegative  plasmas  where  the  attaching 
species  is  consumed  (e.g.,  exdmer  lasers)  the  constriction  described  above  causes  a 
nonuniform  consumption  of  the  attaching  gas  and  may  cause  an  operating  point 
instability  (Ref.  3).  This  effect  alone  can  cause  arcing  and  streamers.  There  may 
also  be  a  measurable  secondary  effect  on  the  self-sustaining  E/N  because  a  plasma 
which  is  constricted  into  a  smaller  region  has  a  higher  resistance. 
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